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1. Introduction 

We cover some useful techniques in computational aspects of ana- 
lytic number theory, with specific emphasis on ideas relevant to the 
evaluation of L-functions. These techniques overlap considerably with 
basic methods from analytic number theory. On the elementary side, 
summation by parts, Euler-Maclaurin summation, and Mobius inver- 
sion play a prominent role. In the slightly less elementary sphere, we 
find tools from analysis, such as Poisson summation, generating func- 
tion methods, Cauchy's residue theorem, asymptotic methods, and the 
fast Fourier transform. We then describe conjectures and experiments 
that connect number theory and random matrix theory. 



2. Basic methods 

2.1. Summation by parts. Summation by parts can be viewed as a 
discrete form of integration by parts. Let / be a function from Z+ to M 
or C, and g a real or complex valued function of a real variable. Then 



(1) Yl /(^)^(^) = E /(^) - / E /(^) 9'm. 

l<n<x \l<n<2: / \l<n<t / 

Here we are assuming that g' exists and is continous on [1, x]. One veri- 
fies this identity by writing the integral as + + . . . + J|^j , noticing 
that the sum in each integral is constant on each open interval, in- 
tegrating, and telescoping. Although our integral begins at t = 1, it 
is sometimes convenient to start earlier, for example at t = 0. This 
doesn't change the value of the integral, the sum in the integrand be- 
ing empty if t < 1 . Formula (QJ) can also be interpreted in terms of the 
Stieltjes integral. 
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A slightly more general form of partial summation is over a set 
{Ai, A2, . . .} of increasing real numbers: 

E f{n)g{Xn) = fin)] g{x) - f [ J] f{n)\ g'{t)dt. 

X„<X \Xn<X J -^-^1 \X„<t J 

As an application, let 

p<x 

denote the number of primes less than or equal to x, and 

p<x 

denote the number primes up to x with each prime weighted by its 
logarithm. The famous equivalence between n{x) ~ x/ log x and 9{x) ~ 
X can be verified using partial summation. Write 

7r(x) = y^logp-^ = 9(x)-^+ f e(t)—^^, 

from which it follows that if 9{x) ~ x then 7r(a;) ~ a;/ log a;. The 
converse follows from 

p<x t 

2.2. Euler-Maclaurin summation. A powerful application of par- 
tial summation occurs when the function /(n) is identically equal to 
1 and the function g{t) is many times differentiable. In that case, 
summation by parts specializes to the Euler Maclaurin formula which 

involves one summation by parts with /(n) = 1 followed by repeated 
integration by parts. For a, 6 e Z,a < b, partial summation gives 

{[t\-a)g\t)dt = bg{b)-ag{a)- / [t\g'{t)dt. 

Here, we have chosen to start the integral at t — a, rather than at 
t — a + 1. Writing [t\ —t — {t}, with {t} the fractional of t we get 

/b pb 
g{t)dt+ / {t}g'{t)dt. 

The second term on the r.h.s. should be viewed as the necessary cor- 
rection that arises from replacing the sum on the left with an integral. 
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The next step is to write {t} = 1/2 + ({t} — 1/2), the latter term 
having nicer properties than {t}, for example being odd and also having 
zero constant term in its Fourier expansion. So 

(2) Yl 9{n) = f g{t)dt+\{g{h)-g{a)) + f\{t} - l/2)g'{t)dt. 

a<n<b "^'^ 

Integrating the second integral repeatedly by parts leads naturally to 
the introduction of Bernoulli polynomials, named after Jacob Bernoulli 
(1654-1705), who discovered them in connection to the problem of 
studying sums of positive integer powers of consecutive integers. Dur- 
ing the 1730's Euler (1707-1783), who studied mathematics from Ja- 
cob's brother Johann (1667-1748), developed the summation formula 
being described in connection with computing reciprocals of powers 
and Euler 's constant. 

2.2.1. Bernoulli Polynomials. The Bernoulli polynomials are defined 
recursively by the following relations 

Bo{t) = 1 

B'.it) = kBk^iit), k>l 

»i 

Bkit)dt = 0, k>l. 



The second equation determines Bk(t) recursively up to the constant 
term, and the third equation fixes the constant. The first few Bernoulli 
polynomials are listed in Table H 



k 


Bk{t) 





1 


1 


t- 1/2 


2 


-t+l/Q 


3 


- 3/2t^ + l/2t 


4 


t^-2t^ + t^- 1/30 


5 


- 5/2t^ + 5/3t^ - l/6t 



Table 1. The first few Bernoulli polynomials 



Let Bk = -Bfc(O) denote the constant term of Bk(t). Bk is called 
the k-th Bernoulli number. We state basic properties of the Bernoulli 
polynomials. Expansion in terms of Bernoulli numbers: 

Bk{t) = J2('')Bk^rnt"', k>0 
^"^^ 



4 MICHAEL RUBINSTEIN 

Generating function: 



- 1 

Fourier series: 



J2Bk{t)z''/k\, \z\<27i 



(3) = t(Z 

1 



(4) = ,>2. 

Functional equation: 

5fc(t) = (-l)^5fc(l-t), A;>0 
Difference equation: 

(5) =t^ A:>o 

Special values: 




k odd, > 3 
= 1 

I.e. 

(6) Bk{l) = Bk{0), unless k = 1 

Recursion: 

k-l 



j:Qb^=o, k>2 



in=0 

Equation Q can be obtained directly. The other formulae can be 
verified using the defining relations and induction. 

Property (jlj) can be used to obtain a formula for ({2m). Let 



1 



Taking t = 0, k = 2m, even, in the Fourier expansion of Bk{{t}) gives 

(-1)-+H2m)! 
^2m ^^^^p;;^ 2C(2m) 

so that 

(-l)™+i(27r)2™ 
^'^'"> = 2(2n;). 
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a formula discovered by Euler. Because C(2^) — >■ 1 as m — > oo, we 
have 

(-l)™+i2(2m)! 



2m 



(27r)2"^ 
as m — >■ oo. 

2.2.2. Euler- Maclaurin continued. Returning to (j2I), we write 
\{t}-l/2)g'{t)dt= t B,{{t})g'{t)dt. 



Breaking up the integral = J^^^ + J^^i + ■ ■ ■ JI^_i, integrating by 
parts, and noting that -82(1) = -82(0), we get, assuming that g^'^^ exists 
and is continous on [a, b], 

Repeating, using -Bjt(l) = Bk{0) if > 2, leads to the Euler-Maclaurin 
summation formula. Let i^' be a positive integer. Assume that g^^^ 
exists and is continous on [a,b]. Then 



9{n)= / g{t)dt + Y, 

a<n<b k=l 



k\ 

BKm)g^''\t)dt. 



4)^^+1 rb 



K\ 

2.2.3. Application: Sums of consecutive powers. We apply Euler-Macluarin 
summation to obtain Bernoulli's formula for sums of powers of consec- 
utive integers. Let r > be an integer. Then 

^ i?,+i(iV + l)-i?,+i(l) 
r + 1 

n=l 

We can verify this directly using property (0), substituting n = 1,2, . . . , N, 
and telescoping. However, it is instructive to apply the Euler-Maclaurin 
formula, which, once begun, carries through in an automatic fashion. 
In this example, we have g{t) = t*". Notice that g^^~^^\t) = 0, and that 



r(r — 1) . . . (r — m + 1)A^'' m < r — 1 
0, m > r. 
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If m = we set r(r — 1) . . . (r — m + 1) = 1. Then 
^n'^ = / fdt + Y, , , V (r - 1) ■ ■ ■ (r - A; + 2)Ar'' 

n=l A:=l 



A;! 



Jo Z-.r-k + l\kJ 

J = J {-^y Br{-t)dt = j Br{t+l)dt 



k=0 

= (5,+i(iV + l)-5,+i(l))/(r + l). 
If r > 1, the last hne simphfies according to © and equals 

Br+ljN + 1) - Br+l 

r + 1 

2.2.4. Application: ({s). The Euler-Maclaurin formula can be used to 
obtain the analytic continuation of ({s) and also provides a useful ex- 
pansion for its numeric evaluation. Consider 



N N 

n 

1 2 

with 3fJs > 1. We have started the sum at = 2 rather than n = 
1 to avoid difficulties near t = below. Applying Euler-Maclaurin 
summation, with g{t) = t~^, g^"^\t) = (— l)'"s(s + 1) . . . (s + m — 
we get 

1 fc=i ^ ^ 



K 



^ j\Ki{t})t-'-''dt. 



Evaluating the first integral, taking the limit as — oo, with 9fJs > 1, 

we get 

(7) 

cw = — +5+E( )f-( K 

While we started with 3?s > 1, the r.h.s. is meromorphic for > 
—K + 1, so gives the meromorphic continuation of ({s) in this region, 
with the only pole being the simple pole at s = 1. 
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I — n N / 



Taking s = 2- K, K >2, 

[ k j''^^- 
Thus, 

^(l_2m) = -52^/(2m), m= 1,2,3,... 
C(-2m) = 0, m = 1,2,3,... 
C(0) = -1/2. 

Applying the functional equation for (see for example Roger Heath- 
Brown's notes) 

7r-/2r(s/2)C(s) = 7r-(i-^)/2r((l - s)/2)C(l - s) 

and 

r(i/2) = .'/'^ = ^^^...d5|^r(i/2-™) 

gives another proof of Euler's identity 

(2m = ^ ' \/ — B2m, m> 1. 

' 2(2m)! 

2.2.5. Computing Q{s) using Euler-Maclaurin summation. Next we de- 
scribe how to adapt the above to obtain a practical method for numeri- 
cally evaluating C,{^)- From a computational perspective, the following 
works better than using ((7j). Let iV be a large positive integer, pro- 
portional in size to We will make this more explicit shortly. For 
> 1, write 

00 N oo 

(8) C(s) = 5^n- = 5^n- + 5^n-. 

1 1 N+l 

The first sum on the r.h.s. is evaluated term by term, while the second 

sum is evaluated using Euler-Maclaurin summation 

(9) 

As before, the r.h.s. above gives the meromorphic continuation of the 
l.h.s. to > —K + 1. Breaking up the sum over n in this fashion 
allows us to throw away the integral on the r.h.s., and obtain sharp 
estimates for its neglected contribution. First, from property (jH), 
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It is convenient to take K = 2Ko, even, in which case we have from 



Therefore, for s 



< 



\B2Kom)\ < 

= a + ir, a > -2Ko + 1, 
S + 2K0- 1 
2Ko 

j^-<T-2Ko+l 



B2K,{{t})t 



-S-2K, 



N 



s + 2Kq 
2Ko 
S + 2K0- 



1 



2Ko 



1 



11 



(T + 2K0-I 

\s + 2Ko-l\ 
a + 2KQ-l 



(y + 2Ko 
S + 2K0- B2K 
2K0-I , 

last term taken |. 



2K^ 



A more precise estimate follows by comparison of -82/^0 with (^(2K^ 
and we have that the remainder is 



< 



n 

i=o 



l-s + j| 

277 



We start to win when 27rA^ is bigger than |s|,|s + l|,...,|s + 2Kq ~ 2 | . 
There are two parameters which we need to choose: Kq and A^, and we 
also need to specify the number of digits accuracy, Digits, we desire. 
For example, with a > 1/2, taking 



with 



27rA^ > 10|s + 2A:o 



27^0 - 1 > Digits + - logiods + 27^0 - 1|) 



achieves the desired accuracy. The main work involves the computation 



of the sum consisting of 0(|s|) terms. Later we will examine 

the Riemann-Siegel formula and its smoothed variants which, for C^i^s)^ 
involves a main sum of 0(|s|^/^) terms. However, for high precision 
evaluation of C^i^s)^ especially with s closer to the real axis, the Euler- 
Maclaurin formula remains an ideal method allowing for sharp and 
rigorous error estimates and reasonable efficiency. 

In fact, we can turn the above scheme into a computation involving 
0(|s|^/^) operations but requiring 0( (Digits + log log precision 
due to cancellation that occurs. In dHl) choose A^ ~ |10s/(27r)|^/^ and 
assume that 3fJs > 1/2. Expand i^^o its Fourier series (jH). 

We only need M = 0{\s\^^'^) terms of the Fourier expansion to as- 
sure a contribution from the neglected terms smaller than the desired 
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precision. Each term contributes 

(10) —^—j. r e2-"^*t-^-^rft, 



N 



K ) (2TiimY 
so the neglected terms contribute altogether less than 

1 \s^K-\\ \s+j 



n 




N'^ a + K-l \ J-l 2nN 

\ j=o 

Here we have combined the ±m terms together. Comparing to an 
integral, the sum above is < 2/{{K — 1)M^~^) and so the neglected 
terms contribute less than 

2 \s + K-l\^ \s + 3\ 



{K - 1)N'' a + K -1 J-^-L 27rMN' 

We start to win when 2itMN exceeds |s|, . . . , |s + — 2|. For a > 1/2, 
choose K > Digits + logio(|s + K - 1\) + 1 and M = N with 

27iMN > 10|s + A'-2|. 

Asymptotically, we can improve the above choices so as to achieve 
M = N ^ |s|i/V(27r), the same as in the Riemann-Siegel formula. The 
only drawback is that extra precision as described above is needed. The 
individual terms summed in (jH)) are somewhat large in comparison to 
the final result, this coming form the binomial coefficients which have 
numerator (s + — 2) . . . (s + l)s, and this leads to cancellation. 

Finally to compute the contribution to the Fourier expansion from 
the terms with |m| < M, we assume that 41^^ so that the terms ±m 
together involve in (jlUj) the integral 

/■oo 



cos{2TTmt)r'~^dt = {2TTmy+^-^ / cos{u)u-'-^ du. 

N J2-KmN 

This can be expressed in terms of the incomplete F function 



oo 
COS 

w 



See Section 3 which describes properties of the incomplete F function 
and methods for its evaluation. 

The Euler-Maclaurin formula can also be used to evaluate Dirichlet 
L-functions. It works in that case due to the periodic nature of the 
corresponding Dirichlet coefficients. For general L-functions, there are 
smoothed Riemann-Siegel type formulae. These are described later. 
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2.3. Mobius inversion with an application to sums and prod- 
ucts over primes. Computations in analytic number theory often in- 
volve evaluating sums or products over primes. For example, let 7r2(x) 
denote the number of twin primes {p,p + 2), with p and p + 2 both 
prime and less than or equal to x. The famous conjecture of Hardy 
and Littlewood predicts that 

TTs X ~ 2 11- — -. 

Generally, it is easier to deal with a sum rather than a product, so we 
turn this product over primes into a sum by expressing it as 

exp I J2 log(l - 2/p) - 2 log(l - 1/p) I . 

\P>2 / 

Letting f{p) = log(l — 2/p) — 21og(l — l/p), we have 



oo 



/(P) = -E " 



m=l 



hence 



(11) f^p) = - E - 1/2^ 

with 



m 

p>2 m=l 



h{s) = J2p''^ > 1- 

p 

We therefore need an efficient method for computing h{m). This will 
be dealt with below. Notice that h{m) — 1/2'" ~ 1/3"^ so the sum on 
the r.h.s. of (jllj) converges exponentially fast. We can achieve faster 
convergence by writing 

p>2 2<p<P p>P 

summing the terms in the first sum, and expressing the second sum as 

-E^;^(M"^)-V2'"-...-i/p'"). 

m=l 

A second example involves the computation of constants that arise in 
conjectures for moments of C{s). The Keating-Snaith conjecture |KeSj 
asserts that 

(12) M,(T) := 1 |C(l/2 + ^t)\"'dt ~ ^(logT)^^ 
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1 \ \ ^ I m + k — 1 



m 



?Ti=0 

n(-^) gcr 



and 



9k = n 



,.0 (k+jy- 

The placement of fc^! is to ensure that g^. is an integer |CFj . Keating 
and Snaith also provide a conjecture for complex values, > —1/2, 
of which the above is a special case. Keating and Snaith used random 
matrix theory to identify the factor gk- The form of (fT^ . without 
identifying gk, was conjectured by Conrey and Ghosh jCGi . 

The above conjecture gives the leading term for the asymptotics for 
the moments of \({l/2 + it)\. In |CFKRS] a conjecture is given for the 
full asymptotics of MkiT): 

M,(T)~5^c.(fc)(logT)'='-^ 

r=0 

where Co(A;) = akgk/k^^- coincides with the Keating-Snaith leading term 
and where the degree k"^ polynomial is given implicitly as an elabo- 
rate multiple residue. Explicit expressions for Cr{k) are worked out 
in jCFKRSS] and are given as co{k) times complicated rational func- 
tions in k, generalized Euler constants, and sums over primes involving 
log(p), 2Fi{k, k, l;p~^) and its derivatives. One method for computing 
the Cj.(/c)'s involves as part of a single step the computation of sums of 
the form 

(14) 5^^^^^, m = 2,3,4,... r = 0,l,2,.... 

p ^ 

We now describe how to efficiently compute h{s) = "YlipP'^ the 
sums in fll4j) . Take the logarithm of 

c(s)=n(i-^^~')"'' ^^>i 

p 

and apply the Taylor series for log(l — x) to get 

(15) logC(s) = V —hims), > 1. 

m 

m=l 
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Let fi{n), the Mobius ii function, denote the Dirichlet coefficients of 
1/C(^): 

oo 



We have 
fi{n) - 
and 



if n is divisible by the square of an integer > 1 

^ number of prime factors of n jf 77, jg SCJUarefree 



1 if r = 1 
otherwise. 



n|r 

The last property can be proven by writing the sum of the left as 
np|r(^ — 1), and it allows us to invert equation (fT3j) 

a(m) , ^, , aim) h(mns) 

m ^-^ m ^-^ n 

m=l m=l n=l 

r=l m\r 



I.e. 



(16) EP"^ = E^l°gC(m.). 

p m=l 

This is an example of Mobius inversion, and expresses h{s) as a sum 
involving (. Mobius inversion can be interpreted as a form of the sieve 
of Eratosthenes. 

Notice that C{i^s) = 1 + 2~™"^ + 3-"*^ + . . . tends to 1, and hence 
log C(ms) tends to 0, exponentially fast as m 00. Therefore, the 
number of terms needed on the r.h.s. of (fT^ is proportional to the 
desired precision. 

To compute the series appearing in (fT^ we can differentiate h{s) 
r times, obtaining 

(17) ^(!2ii^ = (-irf ^^M(,ogc(,„.))M. 

p m=l 

In both ()16|) and (fTTj) . we can use Euler-Maclaurin summation to 
compute C and its derivatives. The paper of Henri Cohen |Cj is a good 
reference for computations involving sums or products of primes. 
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2.4. Poisson summation as a tool for numerical integration. 

Let / e L^{R) and let 



/oo 
-oo 



denote its Fourier transform. The Poisson summation formula asserts, 
for /, / e L^(R) and of bounded variation, that 

oo oo 

E fin) = Yl />)■ 

n=— oo n=— oo 

We often encounter the Poisson summation formula as a potent the- 
oretical tool in analytic number theory. For example, the functional 
equations of the Riemann ^ function and of the Dedekind rj function 
can be derived by exploiting Poisson summation. However, Poisson 
summation is often overlooked in the setting of numerical integration 
where it provides justification for carrying out certain numerical inte- 
grals in a very naive way. 

Let A > 0. By a change of variable 

oo oo 

A J2 /(^A)= J2 /(n/A) = /(0) + J]/(n/A) 

n=— oo n=— oo n^O 



so that 



J — c 



f{t)dt-A Y /(nA) = -5^/(n/A) 

n=— oo n^O 

tells us how closely the Riemann sum A Yl'^=_^ f{nA) approximates 
the integral f{t)dt. 

The main point is that if / is rapidly decreasing then we get enormous 
accuracy from the Riemann sum, even with A not too small. For 
example, with A = 1/10, the first contribution comes from /(±10) 
which can be extremely small if / decreases sufficiently fast. 

As a simple application, let f{t) = exp(— Then f{y) = \/27rexp(— 27r^|/^), 
and so 

5;/(n/A)=0(exp(-27rVA^)). 

Therefore 

/oo °° 
exp(-tV2)dt - A Y exp(-(nA)V2) = 0(exp(-27rVA2)). 
-oo , „ 
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As everyone knows, the integral on the Lh.s. equals v^27r. Taking 
A = 1/10, we therefore get 

oo 

A exp(-(r2A)V2) = v^ + e 



n=— oo 

-857 



with e ~ 10 . We can truncate the sum over n roughly when 



2 ' A2' 

i.e. when n > 271/ A"^. So only 628 terms (combine ±n) are needed to 
evaluate V2tt to about 857 decimal place accuracy! 

This method can be applied to the problem of computing certain 
probability distributions that arise in random matrix theory. Let U be 
an N X N unitary matrix, with eigenvalues exp(i^i), . . .exp{i6^), and 
characteristic equation 

JV 

1 

evaluated on the unit circle at the point exp(z^). In making their 
conjecture for the moments of |C(l/2 + it)\, Keating and Snaith |KeSj 
studied the analogous random matrix theory problem of evaluating the 
moments of \Z{U,6)\, averaged according to Haar measure on U(A^). 
The characteristic function of a matrix is a class function that only 
depends on the eigenvalues of the matrix. For class functions, the Weyl 
integration formula gives Haar measure in terms of the eigenangles, the 
invariant probability measure on U(A^) being 



^ ' ' l<j<m<N 



N- 



Therefore, Mi\f{r), the rth moment of \Z{U,6)\, is given by 

Mjv(r) = — ^ r... r n \e''^-e''-\'\ZiU,9)\^d9,...de^, 

^ ^ ■ "^0 l<j<m<N 

for 3?r > —1. This integral happens to be a special case of Selberg's 
integral, and Keating and Snaith consequently determined that 

Notice that this does not depend on 9. 

Say we are interested in computing the probability distribution of 
\Z(U,9)\. One can recover the probability density function from the 



METHODS AND EXPERIMENTS 



15 



moments as follows. We can express the moments of \Z{U, 6) \ in terms 
of its probability density function. Let 

prob(0 < a < < 6) = f pN{t)dt. 

J a 

Then 

/•oo 

(18) M^(r) = / pN{t)fdt 



Jo 

is a Mellin transform, and taking the inverse Mellin transform we get 

with u to the right of the poles of MAr(r), u > —1. There is an extra 
1/t in front of the integral since the Mellin transform ()18p is evaluated 
at r rather than at r — 1. 

To compute pAr(t) we could shift the line integral to the left picking 
up residues at the poles of Mjv(r), but as grows this becomes bur- 
densome. Instead, we can compute the inverse Mellin transform (jl9j) 
as a simple Riemann sum. 

Changing variables we have 



Let 



f = ^ n r(j)r(j + t/ + 2|/) 

This function also depends on u and N, but we do not include them 
explicitly on the l.h.s. so as to simplify our notation. The above integral 
equals 

/oo 
ft{y)dy. 
-oo 

To estimate the error in computing this integral as a Riemann sum 
using increments of size A, we need bounds on the Fourier transform 

/oo 
/*(2/)e-2™^rf2/. 
-oo 

However, 

/i(l/)e-2™^ = /te2.u(y)e2™('^+i) 

and so 



16 



MICHAEL RUBINSTEIN 



Now, pN(t) is supported in [0,2^], because < \Z{U,9)\ < 2^. Hence 
if M > (Ariog2-logt)/(27r) then ft{u) = 0. Thus, for < t < 2^, if we 
evaluate ()20|) as a Riemann sum with step size A < 27r/(A^log2 — logt) 
the error is 

J2ft{n/A) = J2kn/A) 

since the terms with n > are all zero. On the other hand, with n < 
we get 

where Pmax denotes the maximum of PAr(t) (an upper bound for pAr(t) 
can be obtained from (fT^ ). 
Therefore, choosing 

A = 

Digits log 10 + log 2 - log t 

and setting u = we have 

fti-H/A) < (10-^^s'*^2-^t)l-lp^ax 
Summing over n = —1,-2,— 3,... we get an overall bound of 

We could choose u to be larger, i.e. shift our line integral (fTT?|) to 
the right, and thus achieve more rapid decay of as -u — — oo. 

However, this leads to precision issues. As u increases, the integrand 
in ()19|) increases in size, yet PAr(t) remains constant for given N and t. 
Therefore cancellation must occur when we evaluate the Riemann sum 
and higher precision is needed to capture this cancellation. We leave 
it as an excercise to determine the amount of precision needed for a 
given value of u. 

Another application appears in |RSj where Poisson summation is 
used to compute, on a logarithmic scale, the probability that 7r(a;), the 
number of primes up to x, exceeds Li(a;) = dt/\og(t). The answer 
turns out to be .00000026 . . . 

Later in this paper, we apply this method to computing certain com- 
plicated integrals that arise in the theory of general L-functions. 

3. Analytic aspects of L-function computations 

3.1. Riemann-Siegel formula. The Riemann Siegel formula expresses 
the Riemann ( function as a main sum involving a truncated Dirich- 
let series and correction terms. The formula is often presented with 
3?s = 1/2, but can be given for s off the critical hne. See |USj for a 
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nice presentation of the formula for 1/2 < 3?s < 2 and references. Here 
we stick to 3fts = 1/2. 
Let 

Z{t) = e'^^'\{l/2 + it) 

The rotation factor e*^'-*^ is chosen so that Z{t) is reaL 

For t > 271, let a = (t/(27r))i/^ = [a\, p = {a} = a - [a\ the 
fractional part of a. Then 

N 

Z{t) = 2J2 n'^^^ cos(t log(n) - e{t)) + R{t) 



n=l 



where 



with 



^^+1 — Crip) 



,1/2 

r=0 



(:7o(p) = ^(p) := cos(27r(p2 _ p _ 1/16))/ cos(27rp) 
1 

96^2- 



Clip) = -7^.^''\P) 
1 



C2(p) = — ^—ri^^'^Hp) + — ^i^^^Hp)- 

^'^^ 184327r4^ ^'' Uti^^ ^'^^ 

In general [Ej, Cj{p) can be expressed as a linear combination of the 
derivatives of if). We also have 

Rmit) = 0(t-(2™+3)/4)_ 

Gabcke [G] showed that 

|i?i(t)| < .053r^/^ t> 200. 

The bulk of computational time in evaluating ({s) using the Riemann- 
Siegel formula is spent on the main sum Xl^Li ^"^^^ cos(t log(n) — ^(t)). 
Odlyzko and Schonhage |()Sj [O] developed an algorithm to compute 
the main sum for T <t < T + T^^^ in 0{t^) operations providing that a 
precomputation involving ©(T^/^+e^ operations and bits of storage are 
carried out beforehand. This algorithm lies behind Odlyzko's monu- 
mental ( computations [0] |U2j . An earlier implementation proceeded 
by using the Fast Fourier Transform to compute the main sum and 
its derivatives at equally spaced grid points to then compute the main 
sum in between using Taylor series. This was then improved [01 4.4] 
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to using just the values of the main sum at equally spaced points and 
an interpolation formula from the theory of band-limited functions. 

Riemann used the saddle point method to obtain Cj, for j < 5. 
The reason that a nice formula works using a sharp cutoff, truncating 
the sum over n at A^, is that all the Dirichlet coefficients are equal to 
one. Riemann starts with an expression for ({s) which involves the 
geometric series identity 1/(1 — x) = X]^") Taylor coefficients on 
the right being the Dirichlet coefficients of C{s). For general L-functions 
smoothing works better. 

3.2. Smoothed approximate functional equations. Let 

L(.) = vM 

n=l 

be a Dirichlet series that converges absolutely in a half plane, 3?(s) > 
(Ti, and hence uniformly convergent in any half plane 9fJ(s) > (72 > 0"i 
by comparison with the series for L((j2). 
Let 

(23) K{s)=Q^{^T{k,s + \,)^L{s), 

with Q,K,j G M^, 3?Aj > 0, and assume that: 

(1) A(s) has a meromorphic continuation to all of C with simple 
poles at si, . . . , and corresponding residues ri, . . . , r^. 

(2) (functional equation) A(s) = ujK{1 — s) for some G C, 7^ 0. 

(3) For any a < j3, L{a + it) = O(expt^) for some A > 0, as 
\t\ ^ 00, a < a < (3, with A and the constant in the 'Oh' 
notation depending on a and /3. 

Remarks . a) The 3rd condition, L(cr + it) = 0{expt^), is very 
mild. Using the fact that L{s) is bounded in > a2 > cri, the 
functional equation and the estimate \29^) . and the Phragmen-Lindelof 
Theorem |Rudj we can show that in any vertical strip a < a < (3 , 

L{s) = 0{t^), for some b > 

where both b and the constant in the 'Oh' notation depend on a and (3. 

b) Ifb{n), Xj G M, then the second assumption reads A(s) = ujA{1 — s). 

c) In all known examples the nj 's can be taken to equal 1/2. It is useful 
to know the Legendre duplication formula 

(24) T{s) = {2tt)-^/^2'~^/^T{s/2)T{{s + l)/2). 

However, it is sometimes more convenient to work with \2'J^) . and we 
avoid specializing prematurely to Kj = 1/2. 
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d) The assumption that L{s) have at most simple poles is not crucial 
and is only made to simplify the presentation. 

e) From the point of view of computing A(s) given the Dirichlet co- 
efficients and functional equation, we do not need to assume an Euler 
product for L{s) . Without an Euler product, however, it is unlikely that 
L{s) will satisfy a Riemann Hypothesis. 

To obtain a smoothed approximate functional equation with desir- 
able properties we introduce an auxiliary function. Let g : C C he 
an entire function that, for fixed s, satisfies 

\A{z + s)g{z + s)z''^\ 

as \'^z\ — > oo, in vertical strips, —a <^z<a. The smoothed approx- 
imate functional equation has the following form. 

Theorem 1. For s ^ {si, . . . , Se}, and L{s), g{s) as above, 
Hs)g{s) = 2_^ + Q y^fi{s,n) 

(25) +^Qi-^M/,(l_,,^) 



n=l 



where 



/i(s,n) = — / Y[T{Kj{z + s) + \j)z-^g{s + z){Q/nydz 

(26) 

f,{l-s,n) = — / l[TiK,iz+l-s)+Y,)z-'g{s-z){Q/nydz 

■^'^'^ J u-ioo 

with V > max {0, — 3?(Ai/ki + s), . . . , — 3?(Aa/ + s)}. 

Proof. Let C be the rectangle with verticies {—a, —iT), (a, —iT), (a, iT), 
(—a, iT), let s G C — {si, . . . , s^}, and consider 



(27) 



/ A{z + s)g{z + s)z ^dz. 
Jc 



(integrated counter-clockwise), a and T are chosen big enough so that 
all the poles of the integrand are contained within the rectangle. We 
will also require, soon, that a > ai — 3fJs. On the one hand ^T7\ equals 



(28) A(sMs) + ±rMM 



k=l - ' 
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since the poles of the integrand are included in the set {0, si — s, . . . , — s}, 
and are all simple. Typically, the set of poles will coincide with this 
set. However, if A{s)g{s) = 0, then 2 = is no longer a pole of the inte- 
grand. But then A{s)g{s) contributes nothing to (PHj) and the equality 
remains valid. And if g{sk) = 0, then there is no pole at z = — s but 
also no contribution from rkg{sk)/{sk — s). 

On the other hand, we may break the integral over C into four inte- 
grals: 

= +/ +/ +/ 

JC Ja-iT J a+iT J ~ a+iT J-a-iT 



Ci J C2 J C3 J Ca 



The integral over Ci, assuming that a is big enough to write L{s + z) 
in terms of its Dirichlet series i.e. a > ai — is 



h(n) 1 

^^E^^ / l[^{'^^(' + s) + X,)z-'g{s + z){Q/nrdz. 

n=l ^ j=l 

We are justified in rearranging summation and integration since the 
series for L{z+s) converges uniformly on Ci. Further, by the functional 
equation, the integral over C3 equals 



27n 



/-a~tl 
A(l - 'zT^)g{z + s)z~^dz 
-a+iT 

'Q'"T.V-sV- / ll^{'^^^^-s-z) + X,)z-'g{s + z){Q/nr^dz 

, 1 " J-a+iT „_i 



-a+iT 

°° lif^\ 1 p-a-iT o. 

1 — e ^ 

n=l j=l 

^Q'"^E4^T- / ll^{'^,a-s + z) + X,)z-'g{s-z){Q/nydz. 

n ZTll Ja-iT ■ 1 

n=l 1=1 



Letting T —>■ 00, the integrals over C2 and C4 tend to zero by our 
assumption on the rate of growth of g{s), and we obtain (j2S)). The 
integrals in ()26p are, by Cauchy's Theorem, independent of the choice 
of z/, so long as > max {0, — 3?(Ai/ki + s), . . . , —^{Xa/na + s)}. 

□ 



3.3. Choice of g{z). Formulae of the form (^3)) are well known jL] [F?j . 
Usually, one finds it in the literature with g{s) = 1. For example, for the 
Riemann zeta function this leads to Riemann's formula jEl pg 179] jTi[ 
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pg 22] 

11 °° 1 

7r-^/'T{s/2)C{s) = -- - ^ + TT-/^ 5^ -r(./2, Trn^) 

n=l 

oo _ 

+^(-i)/2^^r((l-.)/2,vrn^) 

n=l 

where r(s,iu) is the incomplete gamma function (see Section ITlj) . 

However, the choice g{s) = 1 is not well suited for computing A(s) 
as \Q{s)\ grows. By Stirling's formula [OH pg 294] 

(29) |r(s)| ~ (27r)i/2|s|--i/2g-|tk/2^ s = a + it 

as |t| — s> oo, and so decreases very quickly as \t\ increases. Hence, 
with g{s) = 1, the l.h.s. of is extremely small for large |t| and 
fixed a. On the other hand, we can show that the terms on the r.h.s., 
though decreasing as ri — > oo, start off relatively large compared to 
the l.h.s.. Hence a tremendous amount of cancellation must occur on 
the r.h.s. and and we would need an unreasonable amount of preci- 
sion. This problem is analogous to what happens if we try to sum 
exp(— a;) = ^(— x)"'/n! in a naive way. If x is positive and large, the 
l.h.s. is exponentially small, yet the terms on the r.h.s. are large before 
they become small and high precision is needed to capture the ensuing 
cancellation. 

One way to control this cancellation is to choose g{s) equal to 
with |5| = 1 and chosen to cancel out most of the exponentially small 
size of the F factors. This idea appears in the work of Lavrik ^ , and 
was also suggested by Lagarias and Odlyzko |LOj who did not imple- 
ment it since it led to complications regarding the computation of ()26|). 
This method was successfully applied in the author's PhD thesis |Ruj 
to compute Dirichlet L-functions and L-functions associated to cusp 
forms and is used extensively in the author's L-function package |Ru3j 
More recently, this approach was used in the computation of Akiyama 
and Tanigawa |A1] to compute several elliptic curve L-functions. 

In fact when there are multiple F factors it is better to choose a 
different S for each F and multiply these together. For a given s let 

tj = ^{KjS + Aj) 

^ in/2, if \tj\ < 2c/{a7f) 

1 c/(a|tj|), if \tj\ > 2c/{an) 

(30) Sj = exp(i sgn(tj)(7r/2 - Oj)). 
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Here c > is a free parameter. Larger c means faster convergence of 
the sums in (I25|) . but also more cancellation and loss of precision. 

Next, we set 

a 

(31) g{z) := Jj57''^'"^^^ = p6-'. 

Because 6j depends on s, the constants 6 and /5 depends on s. We 
can either use a fresh 6 for each new s value, or else modify the above 
choice of tj so as to use the same tj for other nearby s's. The latter is 
prefered if we wish to carry out precomputations that can be recycled 
as we vary s. For simplicity, here we assume that a fresh 6 is chosen as 
above for each new s. 

The choice of g controls the exponentially small size of the F factors. 
Notice that the constant factor (3 = YYj=i appears in every 

term in (j2Sl), and hence can be dropped from g{z) without any effect 
on cancellation or the needed precision. However, to analyze the size 
of the the l.h.s. of and the terms on the r.h.s. this factor is helpful 
and we leave it in for now, but with the understanding that it can be 
omitted. 

To see the effect of the function g{z) on the l.h.s of ()25p we have, 
by ^ and §^ 

\A{so)g{so)\ ~ *-|L(so)| JJ exp (- 7r/2) JJ exp (-c/a) 

\tj\<2c/TT \tj\>2c/n 

> * ■ |L(so)|exp(-c) 

where 

a 

i=i 

We have thus managed to control the exponentially small size of A(s) 
up to a factor of exp(— c) which we can regulate via the choice of c. We 
can also show that this choice of g{z) leads to well balanced terms on 
the r.h.s. of 

3.4. Approximate functional equation in the case of one F- 
factor. We first treat the case a = 1 separately because it is the sim- 
plest, the greatest number of tools have been developed to handle this 
case, and many popular L-functions have a = 1. 
Here we are assuming that 

A(s) = g^F(7s + A)L(s). 



METHODS AND EXPERIMENTS 23 

According to (|!?T|l we should set 

g{s) = 

(we omit the factor /3 as described following (jSH)) with 
and 



ti = $5(75 + A) 

7r/2, if|ti|<2c/7r 
c/\ti\, if \ti\ > 2c/-n 

5i = exp(z sgn(ti)(7r/2 - 6*1)). 



In that case, the function fi{s,n) that appears in Theorem [T] equals 

r-s i-v+ioo 

fi{s, ^) = ^ / r(7(z + + X)z-' {Q/{n6)r dz 



' u—ioo 



2m 

Now 



T{u + 7S + X)u-^ {QI{n5)Y''' du. 



■yiy—ioo 



POD 

(32) T{v + u)u-^= T{v,t)e~^dt, > 0, ^{v + u) > 

Jo 



where 



T{z,w) = / e^^x^^^dx |argw| < tt 

J w 
/■oo 

= w' e-'^'^x'-^dx, ^{w) > 0. 



T{z,w) is known as the incomplete gamma function. By Mellin inver- 
sion 

Similarly 

/2(1 -s,n) = S-r (7(1 -s) + X, {n/{5Q)f''^ 
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We may thus express, when a = 1 and g{s) = 6^"^, (j2SI) as 

k=i ^ 



(33) 
where 



n=l 

oo 

+ '^{Q5)-~^'^Y.'^{n)n^'^G (7(1 - s) + \ {n/{5Q)f^ 

n=l 



00 

wx ^2— 1 , 



(34) G{z,w)=w~'V{z,w) = j e-'^'^x'-^dx, 3ft(u;) > 0. 

Note, from dHU) with a = 1, we have m^'^< > 0, so both {n5/Qf'^ and 
[n/ {5Q))^^'^ have positive 9ft part. 

3.4.1. Examples. 

1) Riemann zeta function, C("5): the necessary background can be 
found in jTIj. Formula for C(s), is 



00 



7r-^/2r^3/2)C(s)r'^ = -i - ^ G (s/2, 7^252) 



s 1 — s 

n=l 



(35) +r^^G((l-s)/2,7rnV52^ 

n=l 

2) Dirichlet L-functions, L{s,x)'- (see jDl chapter 9]). When x is 
primitive and even, x(~l) = 1; "we get 



TT/ 



r(s/2)L(s, x)r ^ = x{n)G (s/2, vrnV/g) 

(x) 



n=l 

00 

n=l 



and when x is primitive and odd, x(~l) = ~1) '^6 get 

s/2 / \ ^/"^ °° 

(^) r(./2 + 1/2)L(., x)r ^ = 5 - ^ x(n)nG {s/2 + 1/2, vrn^^Vg) 



.T^:^^Yx{n)nG ((1 - s)/2 + 1/2, ^7(5^^)) 



n=l 
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Here, r(x) is the Gauss sum 

m=l 

3) Cusp form L-functions: (see |Ug| ). Let f{z) be a cusp form of 
weight k for SL2(Z), /c a positive even integer: 

(1) f{z) is entire on H, the upper half plane. 

(2) fiaz) = (cz + dff{z), a=(^^ e SL2(Z), z E M. 

(3) limt^oo/(^t) =0. 

Assume further that / is a Hecke eigenform, i.e. an eigenfunction of 
the Hecke operators. We may expand / in a Fourier series 

oo 
n=l 

and associate to f{z) the Dirichlet series 

oo 

1 

We normalize / so that ai = 1. This series converges absolutely when 
3?(s) > 1 because, as proven by Deligne jPelj . 

\an\ < ao{n)n^''-'^/\ 

where <7o{n) := J2d\n ^ ~ 0{n'') for any e > 0. 

Lf{s) admits an analytic continuation to all of C and satisfies the 
functional equation 

A^.) := (27r)-T(. + {k - l)/2)Lf{s) = (-l)^/^AXl - 

With our normalization, Oi = 1, the real since they are eigen- 

values of self adjoint operators, the Hecke operators with respect to 
the Petersson inner product (see |Ugt III-12]). Furthermore, the re- 
quired rate of growth on Lf{s), condition 3 on page ITHj follows from 
the modularity of /. 

Hence, in this example, formula ()33|) is 

oo 

(27r)-T(s + {k- l)/2)Lf{s)6-' = (527r)(^"^)/^ ^ a^G {s + {k - l)/2, 27m6) 

n=l 

f y 1 5^ a„G (1 - s + (A; - l)/2, 2nn/5) 

^ ' 71=1 
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4) Twists of cusp forms: Lf{s,x), X primitive, f{z) as in the previous 
example. Lf{s,x) is given by the Dirichlet series 



Lf{s,x) = 



n 



^(fc-l)/2 
1 

Lf{s, x) extends to an entire function and satisfies the functional equa- 
tion 

A/(^,X) := {^yTis + ik-l)/2)Lf{s,x) 

= (-i)'=/M-i)^A;(i-.,x). 

In this example, formula (j33j) is 

±yr{s + {k-l)/2)Lf{s,x)6~^ = 

Y (^nXin)G {s + {k- l)/2, 27rn6/q) 



^ ' n=l 



(_l)fc/2 ^(^) /2,rV'"'^^'>^ 



-X 



/27r\ ^'''^''^ 

[-^) J2 <'nX{n)G {l~s + {k- l)/2, 2nn/{6q)) 



5) Elliptic curve L-functions: (see [Kn,, especially chapters X,XII]). 
Let E be an elliptic curve over Q, which we write in global minimal 
Weierstrass form 

y'^ + cixy + Csy = + C2X^ + C4X + Cq 

where the c/s are integers and the disciminant A is minimal. 
To the elliptic curve E we may associate an Euler product 



(36) Le{s) = \{{l- a,p-'/'-T' 11(1 - + p- 

p\A pfA 



-2s\ 



where, for p f A, ap = p + 1 — ^Ep{Zp), with jj^Ep{2,p) being the 
number of points (x, y) in Zp x Zp on the curve E considered modulo 
p, together with the point at infinity. When p| A, ap is either 1, —1, or 
0. If p f A, a theorem of Hasse states that \ap\ < 2p^l'^ . Hence, (jHBj) 
converges when 3?(s) > 1, and for these values of s we may expand 
Le{s) in an absolutely convergent Dirichlet series 



00 

a 



(37) L,(s) = J2^n 
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The Hasse-Weil conjecture asserts that Le{s) extends to an entire 
function and has the functional equation 

Ms) := f ^ j Vis + 1/2)Le{s) = -e^E{l - s). 

where is the conductor of and e, which depends on ii^, is either 
±1. The Hasse-Weil conjecture and also the required rate of growth 
on Le{s) follows from the Shimura-Taniyama-Weil conjecture, which 
has been proven by Wiles and Taylor |TWj |Wij for elliptic curves 
with square free conductor and has been extended, by Breuil, Conrad, 
Diamond and Taylor to all elliptic curves over Q [BC DTj . 
Hence we have 



r(. + i/2)L^(.)r^ =\wrA E + ^2, iT^nblN^I-") 

' ^ ' n=l 

(9 \ °o 
5^a„G(l-. + l/2,2W(5iVV2))- 
^ n=l 

6) Twists of elliptic curve L-functions: Le{s,x), x ^ primitive 
character of conductor q, (g, A^) = 1. Here Le{s, x) is given by the 
Dirichlet series 



oo 



nV2 



The Weil conjecture asserts, here, that Le{s) extends to an entire func- 
tion and satisfies 

Ae(s,x) := (^) T{s+1/2)Le{s,x) = -ex(-iV)^A£;(l-s, x)- 

Here N and e are the same as for E. In this example the conjectured 
formula is 

/ ]\t1/2 \ / 9 A \ °° 

j r(. + 1/2)Le{s)5-' = l^-^j g a^xin)G {s + 1/2, 2nnS / (qN^/')) 

( \ / <-) \ 1/2 °o 

-ix(-iV)^ j 5: ^rOi{n)G {1-S + 1/2, 2™/(5giVV2)) . 



We have reduced in the case a = 1 the computation of A(s) to one 
of evaluating two sums of incomplete gamma functions. The r(7s + 
X)5~^ factor on the left of (jH^ and elsewhere is easily evaluated using 
several terms of Stirling's asymptotic formula and also the recurrence 
T{z + 1) = zT{z) apphed a few times. The second step is needed for 
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small z. Some care needs to be taken to absorb the e~'^\'^'<^^^^^\l'^ factor 
of r(7s + A) into the e^\'^'<^^+^)\l'^ factor of . Otherwise our effort to 
control the size of r(7s + A) will have been in vain, and lack of precision 
will wreak havoc. 

To see how many terms in are needed we can use the rough 
bound 

/•oo s)[{(ui) 



l{w)-^{z) + 1' 

valid for 3?(w) > ^{z) — 1 > 0. We have put t = a; — 1 in and have 
used t + 1 < e*. Also, for ^{w) > and ^{z) < 1, 

\G{z,w)\ < ——. 

These inequalities tells us that the terms in decrease exponentially 
fast once n is sufficiently large. 

For example, in equation (jHSj) for ({s) we get exponential drop off 
roughly when 

^7m^6^ » 1. 

But 

^7rn^5^ = Trn^m^ ~ 27m^c/t 
so the number of terms needed is roughly 

» (t/c)i/2_ 

3.4.2. Computing r{z,w) . Recall the definitions 

rC^:,^;) = / e'H'^'^dt, |argiy|<7r 



G{z,w) = w T{z,w). 

Let 

PW 

'y{z,w) ■.= T{z) -T{z,w) = / e~'^x''~^dx, > 0, |argw| < vr 

Jo 

be the complimentary incomplete gamma function, and set 
(38) g{z,w)=w~'-f{z,w)= [ e-'"H'-^dt 



so that G{z,w) + g{z,w) = w^^T{z). The function g{z,w)/T{z) is 
entire in z and w. 

The incomplete F function undergoes a transition when |w| is close 
to \z\. This will be described using Temme's uniform asymptotics for 
T{z,w). The transition explains the difficulty in computing T{z,w) 
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without resorting to several different expressions or using uniform asymp- 
totics. 

A combination of series, asymptotics, and continued fractions are 
useful when \z\ is somewhat bigger than or smaller than \w\. When 
the two parameters are close in size to one another, we can employ 
Temme's more involved uniform asymptotics. We can also apply the 
Poisson summation method described in Section 2, or an expansion 
due to Nielsen. Below we look at a few useful approaches. 

Integrating by parts we get 



g{z,w) = e-^Y. 



W-' 



where 




[z + 3-l) ifj>0; 

ifj=0. 



(The case j = occurs below in an expression for G{z, w)). While this 
series converges for ^ 7^ 0, — 1, — 2, . . . and all w, it is well suited, say if 
3^2; > and |w| < a\z\ with < a < 1. Otherwise, not only does the 
series take too long to converge, but precision issues arise. 
The following continued fraction converges for > 

g{z,w) = 

zw 



w 

Z + 1 + 



(z + l)w 

z + 2 ^ 

2w 

Z + 3 + 



iz + 2)w 

z + A- 



Z + 5 + --- 

The paper of Akiyama and Tanigawa |AT| contains an analysis of the 
truncation error for this continued fraction, as well as the continued 
fraction in (jHIIjl below, and show that the above is most useful when 
\w\ < \z\, with poorer performance as |w| approaches \z\. 
Another series, useful when \w\ << 1, is 



g{z,w) = 



This is obtained from (jHHj) by expanding e~"'* in a Taylor series and 
integrating termwise. As |w| grows, cancellation and precision become 
an issue in the same way it does for the sum e""' = '^{—wy/j^-- 
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Next, integrate G{z,w) by parts to obtain the asymptotic series 
G{z, w) = V -f- + tMiz, w) 



with 



eui^z, w) = ^ — w-G"!^ - M, w). 

— 



This asymptotic expansion works well if \w\ > f3\z\ with (3 > 1 and \z\ 
large. In that region the following continued fraction also works well 

(39) G{z,w) = 



1 - z 

w H 



1 + 



w + 



w + 



1 + 



w + ■ ■ ■ 

Temme's uniform asymptotics for T{z, w) provide a powerful tool for 
computing the function in its transition zone and elsewhere. Following 
the notation in [Tj, let 

Qiz,w) = Viz, w) /Viz) 
\ = w/z 
rf/2 = A- 1 -log A 

where the sign of 1] is chosen to be positive for A > 1. Then 
Q{z,w) = ^eTfc{r]{z/2y/^) + RM 

where 



erfc = I e ^ dt. 



oo 

2 



and Rz is given by the asymptotic series, as ^ — *• oo, 

Cn{ri) 



V/2 °° 



^ ' n=0 
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Here 



ciiv) 



A- 1 T] 

111 1 

^ ~ (A -1)3 ~ (A - 1)2 ~ 12(A- 



?7Cn(?7) = — c„_i(?7) + -7„, n > 1 

dr] A — 1 



with 



■1)"7„ 





being the asymptotic expansion of 

r{z) = {z/{2n)y/'{e/zyT{z). 

The first few terms are 70 = 1, 71 = —1/12, 72 = 1/288, 73 = 
139/51840. The singularities at r/ = 0, i.e. X = 1, z = w, are remov- 
able. Unfortunately, explicit estimates for the remainder in truncating 
()40|) when the parameters are complex have not been worked out, but 
in practice the expansion seems to work very well. 

To handle the intermediate region \z\ ~ \w\ we could also use the 
following expansion of Nielsen to step through the troublesome region 
(41) 

7(z, w + d)= ^{z, w) + w^~'e-^ J2 T^^^ " e-^e, (rf)), \d\ < \w\ 

j=o ^ ^' 

where 



m=0 

A proof can be found in |KM()^I] . This expansion is very well suited, 
for example, for L-functions associated to modular forms, since in that 
case we increment w in equal steps from term to term in ()33p and 
precomputations can be arranged to recycle data. Numerically, this 
expansion is unstable if \d\ is big. This can be overcome by taking many 
smaller steps, but this then makes Nielsen's expansion an inefficient 
choice for ({s) or Dirichlet L-functions. 

In computing (j4H) some care needs to be taken to avoid numeri- 
cal pitfalls. One pitfall is that, as j grows, e~'^ej{d) ^ 1. So once 
|l — e~'^ej{d)\ < 10"°'^'*'^, the error in computation of 1 — e~'^ej{d) is 
bigger than its value, and this gets magnified when we multiply by 
(1 — z)j/{—wy. So in computing ((1 — z)j/{—wy) (1 — e^'^ej{dy one 
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must avoid the temptation to view this as a product of (1 — z)j/ {—wY 
and 1 — e~'^ej{d). Instead, we let 

(1 - z)- 

aj{z,w,d) = ^j-^{l-e-^ej{d)). 
Now, 1 — e~'^ej{d) = e~'^{e'^ — ej{d)), and we get 

aj+i{z,w,d) = aj{z,w,d)^^-^^-^ ^ ^^rf"/m! 

= a,iz, w, d)^^-^^^ (1 - l/Pj{d)) , J = 1, 2, 3, . 

where 

oo 

(3,{d) = J2d"'/ij + ^)m- 

Furthermore 

f3j{d)-l^d/{j + 2), as\d\/j^O. 
Hence, for \w\ ~ \z\, we approximately have (as \d\ /j 0) 



w 



< 1 1 + 1±1] Ml ^ Ml I Ml 



Thus, because \d/w\ < 1, we have, for j big enough, that the above is 
< 1, and so the sum in (PT|) converges geometrically fast, and hence 
only a handful of terms are required. 

One might be tempted to compute the Pj{dys using the recursion 

P,+,id) = iPM)-l)ij + 2)/d 

but this leads to numerical instability. The (3j{d)^s are all equal to 
1 + Od{l/{j + 2)) and are thus all roughly of comparable size. Hence, 
a small error due to roundoff in f3j{d) is turned into a much larger 
error in /3j+i(c?), (j + 2)/ \d\ times larger, and this quickly destroys the 
numerics. 

There seems to be some potential in an asymptotic expression due 
to Ramanujan [HI pg 193, entry 6] 

M 

G{z, w) ~ w-'r{z)/2 + e^"" ^pfc(u; -z + as |z| oo, 

fc=0 

for |w — 2;| relatively small, where Pk{v) is a polynomial in v of degree 
2k + 1, though this potential has not been investigated substantially. 
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We list the first few Pk{v)'s here: 
Po{v) =-v + 2/3 



y3 ^2 



3 3 135 

, , 2v'^ Av 8 

P2[v) = 1 1 

^ 15 9 135 135 2835 

v'^ 7v^ 8v^ 16 f 2 16 V 16 

pJv) = \ \ \ \ 

' 105 45 45 135 405 567 2835 8505 

f9 2v^ 2v^ 8v^ 11 62 32 16 16 v 8992 
~~945~315~315^405^ 405 " 2835 ~ 1215 ^ 1701 ^ 2835 ~ 12629925 

t;" 2v^ 2v^ 43 1;^ 41 1;^ 968 1;^ 68 v* 368 v'^ 

~ ~ 10395 ^ 945 ~ 567 ~ 2835 ^ 2835 ^ 2835 " 42525 ~ 2835 ^ 25515 

138064 f 2 35968 w 334144 

^ 12629925 ~ 12629925 ^ 492567075 
It is wortli noting tliat wlien many evaluations of A(s) are required, 
we can reduce through precomputations the bulk of the work to that 
of computing a main sum. This comes from the identity 

G{z, w) = w~^T{z) — g{z, w). 

The above discussion indicates that, in (jHHj) . we should use g{z, w) and 
this identity to compute G{z,w) roughly when \w\ is smaller than \z\. 
For example, with ({1/2 + it), the region \w\ < \z\ corresponds in 
to \nn^6^\ < |l/4 + it/2| and \TTn^/6^\ < \l/A-it/2\. Because \5\ = 1 
this leads to a main sum consisting of approximately |t/(27r)|^/^ terms, 
the same as in the Riemann-Siegel formula. 

3.5. The approximate functional equation when there is more 
than one F-factor, and = 1/2. In this case, the function fi{s,n) 
that appears in Theorem ^ is 

(42) f,(s,n) = — / llTiiz + s)/2 + X,)z-'iQ/i5n)rdz. 

This is a special case of the Meijer G function and we develop some of 
its properties. 

Let M {(pit); z) denote the Mellin transform of 

POO 

Jo 

We will express nj=ir((-z + s)/2 + Xj)z~^ as a Mellin transform 
analogous to 
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Letting 0i * 02 denote the convolution of two functions 

dt 
T 

we have (under certain conditions on 4>i,4>2) 

M {(f)i * (f)2; z) = M {(j),; z) ■ M {(j)2; z) . 

Thus 

°' POO 

(43) T\M{<j),;z)= (0i*...*0„)(t)t^-idt, 

j=i 

with 

poo roo 

(01*- ■ ■*0a)(^^) = / ••• / Mv/tl)Mil/i2) ■ ■ ■<Pa-l{ta-2/ta-l)<Paita-l)-^ 

Jo Jo 



Now 



n r((^+s)/2+A,-)^-' = n + ^)/2 + A,) (r((z + s)/2 + Xa)z-') . 
j=i \j=i J 

But 

T{{z + s)/2 + A) = M (2e-''f^+'; 2) , 

and p2|l gives 

r((z + s)/2 + A)^-^ = M (r(s/2 + A, t'); z) . 

So letting 

'2e-'^t^^^+' j = l,...a-l; 
r(s/2 + A„t2) j=a, 

and applying Mellin inversion, we find that ()42|) equals 

(44) /i(s,n) = r^(0i*---*0,)(n5/Q), 

where 



(01 * • ■ ■ * 0a)(t;) = / . . . / 2»-i ''^'"''^^"^'^^'^^'"^'^ 

\Jl J h ta-l 

Substituting uj = ^l^^ t| and rearranging order of integration this 



becomes 



/oo 
Ex {xv"^) x'/^+M-Mx, 
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where 

(45) = ^ E 



> 1 ' 
a ^ 
i=\ 



(46) 



a— 1 



So, returning to (011), we find that 

PCO 

h (s, n) = {n5/Qf^ {n/Qf / (x {n5/Qf) x'/^+^-^dx. 



Note that because ()42|1 is symmetric in the Aj's, so is 
Similarly 

/oo 
(n/(5Q))') x(i-^)/2+7i-i^^_ 

Hence, 

Q^f[T{s/2 + X,)L{s)6-^ = Y.'-^ 

■ 1 1 ^ ^ 

3=1 k=l 

oo 

+ Kny^Gx {s/2 + /i, {n6/Qf) 

n=l 

oo 

+ ^m~'~'Y^iny-^G,, ((1 - s)/2 +71, (n/(<5g))^) 



(47) 
with 

/oo 
-Ea (xw) x^~^dx 

{fi and -Ea are given by (jlHI)). 
3.5.1. Examples. When a = 2 

Jo ^ 
(48) = 2ir,,_,, (2(^x)V2) = 2ir,,_,, (2(^x)i/2) , 

K being the i^'-Bessel function, so that Gx is an incomplete integral of 
the i^'-Bessel function. 

Note further that if Ai = A/2, Aa = (A + l)/2 then (gHD is 

2Ki/2 {2{wxY/') = {n'/y{wxf/^) e-2(--)^^^ 



36 MICHAEL RUBINSTEIN 

(see !EM()Tp . so G'(a/2,(a+i)/2) {z, w) = 2(27r) V2(4u;)-T(2z-l/2, 2w'/'), 
i.e. the incomplete gamma function. This is what we expect since, us- 
ing (j21I), we can write the gamma factor r((s + A)/2)r((s + A + l)/2) 
in terms of r(s + A), for which the a = 1 expansion, (j33|) . apphes. 

Maass cusp form L-functions: (background material can be found 
in |Buj ) . Let / be a Maass cusp form with eigenvalue A = 1/4 — v^, 
i.e. A/ = A/, where A = —y'^{d/dx^ + d/dy'^), and Fourier expansion 

/(^) = ^a„yVX(27r|n|y)e2™^ 
with a-n = CLn for all n, or a_„ = —an for all n. Let 

oo 
ri=l 

(absolute convergence in this half plane can be proven via the Rankin- 
Selberg method), and let £ = or 1 according to whether a_„ = a„ or 
a^n = —CLn- We have that 

Kf{s) := 7r-T((s + e + v)/2)T{{s + 6- v)/2)Lf{s) 

extends to an entire function and satisfies 

kf{s) = {-lYkf{l-s). 

Hence, formula (jlZj), for Lf{s), is 

7r-T((s + e + v)/2)T{{s + e- v)/2)Lf{s)6-' = 

oo 

{s/2 + e/2, {n67rf) 

n=\ 

+ ^-^{r^/Sy Yl «n^'G'((.+^)A(s-TT)/2) ((1 - s)/2 + s/2, (nn/df) 

n=l 

where, by (gH)), 

/oo 
K,{2n5TTt)f^'-'dt 

/oo 
K,,{2n7ct/S)t-'^'dt. 
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Next, let 

POO 

{z, w) = w'Gx (z, w)= (t) t'-'dt, 

J w 

POO 

(49) r;,(^) = / E^it)t'-'dt, 

Jo 

pw 

-ix{z,w) = / Ex{t)e-^dt, 
Jo 

with Ex given by ()46j) . 
Lemma 1. 

a 

Tx{z) = l[T{z-^i + X,) 

where ji = ^ EJ=i -^i- 

Proof. Let ipjit) = e~^t^^ , j = 1, . . . , a, and consider 



poo poo -L f ^1 



oo 



a-l 



t 



'o Jo ^1 



(we have put tj = ^/'^Xj). Thus, from 

Ex {v) = V~^{tpi *■■■* 1pa){v), 

and hence (jl^ equals 

which, by (gSl) is nj=i r(z - /i + A^). 







□ 



Inverting, we get 

riy+ioo 



1 PU+'lOO 

Ex{t) = — / Tx{z)t-^dz 



with 1/ to the right of the poles of Fa {z) . Shifting the line integral 
to the left, we can express Ex if) as a sum of residues, and hence 
obtain through termwise integration a series expansion for 7a(;z,w). 
An algorithm for doing so is detailed in |Doj . though with different 
notation. Such an expansion is useful for \w\ « 1. That paper also 
describes how to obtain an asymptotic expansion for Ex it) and hence. 
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by termwise integration, for Fx {z, w), useful for \w\ large in comparison 
to \z\. The paper has, implictly, g{z) = 1 and does not control for 
cancellation. Consequently, it does not provide a means to compute 
L-functions away from the real axis other than increasing precision. 

If one wishes to use the methods of this paper to control for can- 
cellation, then one will have w varying over a wide range of values for 
which the series expansion in |Doj is not adequate. We thus need an 
alternative method to compute Gx {z, w) especially in the transition 
zone \z\ ^ \w\. It would be useful to have Temme's uniform asymp- 
totics generalized to handle Gx {z,w). Alternatively, we can apply the 
naive but powerful Riemann sum technique described in section 2. 

3.6. The functions fi{s,n), f 2(1 — s,n) as Riemann sums. Substi- 
tuting z = V + iu into ^I^i we have 

/i(s,n) = — / T[T{K,{s + v + zu) + Xj)^^^-^^-^{Q/ny+'''du. 

^71" J -00 jJl V + IU 

Let 

h{u) = ^\\ T{k^{s + v + iu) + A,)^i^±^i±^(g/n)-+». 

/TT V + lU 

With the choice of g{z) as in ()3ip . an analysis similar to that follow- 
ing (j7T|) shows that hiy) decays exponentially fast as ?/ —00, and 
doubly exponentially fast as y —>■ 00. Hence, we can successfully eval- 
uate /i(s, n), and similarly /2(1 — s,n) as simple Riemann sums, with 
step size inversely proportional to the number of digits of precision 
required. 

The Riemann sum approach gives us tremendous flexibility. We are 
no longer bound in our choice of g{z) to functions for which ()26|) has 
nice series or asymptotic expansions. For example, we can, with A > 0, 
set 

a 

g{z) = exp{A{z - s)')l[6p\ 

The extra factor exp{A{z — s)^) is chosen so as to cut down on the 
domain of integration. Recall that in n) and /2(1 — s, n), g appears 
as g{s ± {v + iu)), hence exp{A{z — s)^) decays in the integral like 
exp(— Am^). Ideally, we would hke to have A large. However, this would 
cause the Fourier transform h{y) to decay too slowly. The Fourier 
transform of a product is a convolution of Fourier transforms, and the 
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Fourier transform of exp{A{v + iuY) equals 

{■n/Af/"^ exp(7ry(2At; - ny) /A). 

A large value of A leads to to a small 1/ A and this results in poor per- 
formance of h{y). We also need to specify f , for the line of integration. 
Larger v means more rapid decay of h{y) but more cancellation in the 
Riemann sum and hence loss of precision. 

Another advantage to the Riemann sum approach is that we can 
rearrange sums, putting the Riemann sum on the outside and the sum 
over n on the inside. Both sums are finite since we truncate them 
once the tails are within the desired precision. This then expresses, 
to within an error that we can control by our choice of stepsize and 
truncation, A(s) as a sum of finite Dirichlet series evaluated at equally 
spaced points and hence gives a sort of interpolation formula for A(s). 
Details related to this approach will appear in a future paper. 

3.7. Looking for zeros. To look for zeros of an L-function, we can 
rotate it so that it is real on the critical line, for example working with 
Z{t), see (1221), rather than C(l/2 + it). 

We can then advance in small steps, say one quarter the average 
gap size between consecutive zeros, looking for sign changes of this real 
valued function, zooming in each time a sign change occurs. Along 
the way, we need to determine if any zeros have been missed, and, if 
so, go back and look for them, using more refined step sizes. We can 
also use more sophisticated interpolation techniques to make the search 
for zeros more efficient [U]. If this search fails to turn up the missing 
zeros, then presumably a bug has crept into one's code, or else one 
should look for zeros of the L-function nearby but off the critical line 
in violation of the Riemann hypothesis. 

To check for missing zeros, we could use the argument principle and 
numerically integrate the logarithmic derivative of the L-function along 
a rectangle, rounding to the closest integer. However, this is inefficient 
and difficult to make numerically rigorous. 

It is better to use a test devised by Alan Turing |l'uj for C(s) but 
which seems to work well in general. Let N{T) denote the number of 
zeros of ({s) in the critical strip above the real axis and up to height 
T: 

N{T) = \{p = /5 + z7|C(p) = 0,0 < /3 < 1,0 < 7 < T}| . 
A theorem of von Mangoldt states that 

(50) N{T) = ^ log(T/(2vr)) -^ + 1 + S{T) + 0{T-') 
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with 

S{T) = 0{\ogT). 

However, a stronger inequality due to Littlewood and with exphcit 
constants due to Turing |Tuj |Le] is given by 

S{t)dt < 2.3 + .1281og(t2/vr) 



(51) 



for all t2 > ti > IGSvr, i.e. S{T) is on average. Therefore, if we 
miss one sign change (at least two zeros), we'll quickly detect the fact. 
To illustrate this. Table El contains a list of the imaginary parts of the 
zeros of ({s) found naively by searching for sign changes of Z{t) taking 
step sizes equal to two. We notice that near the ninth zero on our list 
a missing pair is detected, and similary near the twenty fifth zero. A 
more refined search reveals the pairs of zeros with imaginary parts equal 
to 48.0051508812, 49.7738324777, and 94.6513440405, 95.8706342282 
respectively. 

It would be useful to have a general form of the explicit inequal- 
ity (j51|) worked out for any L-function. The papers of Rumely |Rumj 
and Tollis |lbj generalize this inequality to Dirichlet L-functions and 
Dedekind zeta functions respectively. 

The main term, analogous to (j50j) . for a general L-function is easy to 
derive. Let L{s) be an L-function with functional equation as described 
in ()23|). Let Nl{T) denote the number of zeros of L(s) lying within 
the rectangle I^SsI < T, < 3?s < 1. Notice here we are considering 
zeros lying both above and below the real axis since the zeros of L(s) 
will not be located symmetrically about the real axis if its Dirichlet 
coefficients b{n) are non-real. 

Assume for simplicity that L(s) is entire. The arguement principle 
and the functional equation for L(s) suggests a main term for Nl{T) 
equal to 

2r, lA^/, /r((l/2 + 2T)fi:, + A,)\\ 

If we assume further that the A^'s are all real, then the above is, by 
Stirling's formula, asymptotically equal to 

Nl{T) ~ — log(g) + V \og{TK,/e) + {k,/2 + A, - 1/2)) . 

j=l V / 

A slight modification of the above is needed if L(s) has poles, as in the 
case of See Davenport jDj chapters 15,16] where rigorous proofs 
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are presented for ((s) and Dirichlet L- functions (the original proof 
due to von Mangoldt). 



j 




N{{t,+t,_i)/2)-j + l 


1 


14.1347251417 


-0.11752 


2 


21.0220396388 


-0.04445 


3 


25.0108575801 


-0.03216 


4 


30.4248761259 


0.01102 


5 


32.9350615877 


-0.01000 


6 


37.5861781588 


-0.05699 


7 


40.9187190121 


0.07354 


8 


43.3270732809 


-0.07314 


9 


52.9703214777 


0.81717 


10 


56.4462476971 


2.01126 


11 


59.3470440026 


2.12394 


12 


60.8317785246 


1.90550 


13 


65.1125440481 


1.95229 


14 


67.0798105295 


2.11039 


15 


69.5464017112 


1.94654 


16 


72.0671576745 


1.90075 


17 


75.7046906991 


2.09822 


18 


77.1448400689 


2.10097 


19 


79.3373750202 


1.82662 


20 


82.9103808541 


1.99205 


21 


84.7354929805 


2.09800 


22 


87.4252746131 


2.03363 


23 


88.8091112076 


1.88592 


24 


92.4918992706 


1.95640 


25 


98.8311942182 


3.10677 


26 


101.3178510057 


4.03517 


27 


103.7255380405 


4.11799 



Table 2. Checking for missing zeros. The second col- 
umn lists the imaginary parts of the zeros of found 
by looking for sign changes of ^(t), advancing in step 
sizes equal to two. The third column compares the num- 
ber of zeros found to the main term of N(T), namely to 
N{T) := (T/(27r))log(T/(27re)) + 7/8, evaluated at the 
midpoint between consecutive zeros, with to taken to be 
0. This detects a pair of missing zeros near the ninth and 
twenty fifth zeros on our list. 
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4. Experiments involving L-functions 

Here we describe some of the experiments that reflect the random 
matrix theory philosophy, namely that the zeros and values of L- 
functions behave like the zeros and values of characteristic functions 
from the classical compact groups |KS2j . Consequently, we are inter- 
ested in questions concerning the distribution of zeros, horizontal and 
vertical, and the value distribution of L-functions. 

4.1. Horizontal distribution of the zeros. Riemann himself com- 
puted the first few zeros of and detailed numerical studies were 
initiated almost as soon as computers were invented. See Edwards [E] 
for a historical survey of these computations. To date, the most im- 
pressive computations for C{s) have been those of Odlyzko [0] |U2j and 
Wedeniwski ^Wj. The latter adapted code of van de Lune, te Riele, and 
Winter |LRWj for grid computing over the internet. Several thousand 
computers have been used to verify that the first 8.5- 10^^ nontrivial ze- 
ros of C,{s) fall on the critical line. Odlyzko's computations have been 
more concerned with examining the distribution of the spacings be- 
tween neighbouring zeros, although the Riemann Hypothesis has also 
been checked for the intervals examined. In [O], Odlyzko computed 175 
million consecutive zeros of C,{s) lying near the 10^°th zero, and more 
recently, billions of zeros in a higher region |U2j . The Riemann- Siegel 
formula has been at the heart of these computations. Odlyzko also 
uses EFT and interpolation algorithms to allow for many evaluations 
of C,{s) at almost the same cost of a single evaluation. 

Dirichlet L-functions were not computed on machines until 1961 
when Davies and Haselgrove |DHj looked at several L(s, x) with con- 
ductor < 163. Rumely Rumj, using summation by parts, computed the 
first several thousand zeros for many Dirichlet L-functions with small 
moduli. He both verified RH and looked at statistics of neighbouring 
zeros. 

Yoshida |Y2j has also used summation by parts, though in a 
different manner, to compute the first few zeros of certain higher degree, 
with two or more F-factors in the functional equation, L-functions. 

Lagarias and Odlyzko |LOj have computed the low lying zeros of 
several Artin L-functions using expansions involving the incomplete 
gamma function. They noted that one could compute higher up in 
the critical strip by introducing the parameter 5, as explained in sec- 
tion ESI but did not implement it since it led to difficulties concerning 
the computation of G{z,w) with both z and w complex. 

Other computations of L-functions include those of Berry and Keat- 
ing |BKj and Paris £j (C("S)), ToUis |Toj (Dedekind zeta functions). 
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Keiper jKej and Spira |Sp| ( Ramanujan r L-function), Fermigier jFj 
and Akivama-Tanigawa |Al'j (elliptic curve L- functions) , Strombergs- 
son jStj and Farmer-Kranec-Lemurell |FKLj (Maass waveform L-functions) , 
and Dokchister |Doj (general L-functions near the critical line). 

The author has verified the Riemann hypothesis for various L-functions. 
These computations use the methods described in sectionEland are not 
rigorous in the sense that no attempt is made to obtain explicit bounds 
for truncation errors on some of the asymptotic expansions and contin- 
ued fractions used, and no interval arithmetic to bound round off errors 
is carried out. Tables of the zeros mentioned may be obtained from the 
author's homepage |Ru4j . These include the first tens of millions zeros 
of all L{s, x) with the conductor of x less than 20, the first 300000 ze- 
ros of Lt[s), the Ramanujan r L-function, the first 100000 zeros of the 
L-functions associated to elliptic curves of conductors 11, 14, 15, 17, 19, 
the first 1000 zeros for elliptic curves of conductors less than 1000, the 
first 100 zeros of elliptic curves with conductor less than 8000, and 
hundreds/millions of zeros of many other L-functions. 

In all these computations, no violations of the Riemann hypothesis 
have been found. 



4.2. Vertical distribution: correlations and spacings distribu- 
tions. The random matrix philosophy predicts that various statistics 
of the zeros of L-functions will mimic the same statistics for the eigen- 
values of matrices in the classical compact groups. 

Montgomery |Moj achieved the first result connecting zeros of ({s) 
with eigenvalues of unitary matrices. Write a typical non-trivial zero 
of C as 

1/2 + Z7. 

Assume the Riemann Hypothesis, so that the 7's are real. Because the 
zeros of ({s) come in conjugate pairs, we can restrict our attention to 
those lying above the real axis and order them 

< 7i < 72 < 73 . . . 

We can then ask how the spacings between consecutive zeros, 7^+1 —ji, 
are distributed, but first, we need to 'unfold' the zeros to compensate 
for the fact that the zeros on average become closer as one goes higher 
in the critical strip. We set 

(52) 7.^7.i25^^« 
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and investigate questions involving the 7's. This normahzation is cho- 
sen so that the mean spacing between consecutive 7's equals one. Sum- 
ming the consecutive differecnces, we get a telespcoping sum 

> (7m - 7*) = l{T) + 0(1) = 7(T) + 0(1) 

^ — ' ZTT 

where 7(T) is the largest 7 less than or equal to T. By dHOI), the r.h.s 
above equals 

iV(7(T)) + 0(log(7(T))) = N{T) + 0(log(T)), 

hence 7j+i — 7^ has mean spacing equal to one. 

From a theoretical point of view, studying the consecutive spacings 
distribution is difficult since this assumes the ability to sort the zeros. 
The tool that is used for studying spacings questions about the zeros, 
namely the explicit formula, involves a sum over all zeros of C{s)i and 
it is easier to consider the pair correlation, a statistic incorporating 
differences between all pairs of zeros. Montgomery conjectured that 
for < a < /? and M 00, 

M-i|{l<2<j<M:7,-7, G[a,/3)}| 

' ' ' sm Trt ^ 



(53) ~ / I 1 - ) I dt. 

Notice that M~^, and not, say, (^2^) ^, is the correct normalization. 
For any j there, are just a handful of i's with 7^ — 7^ G [a, /?). 
Montgomery was able to prove that 

(54) M- E^/(7.-70-/°°/(t)(l-(^)')cit. 

as M — 00, for test functions / satisfying the stringent restriction that 
/ be supported in (—1, 1). 

An equivalent way to state the conjecture as M — > 00, and one which 
Odlyzko uses in his numerical experiments, is to let 

(55) = 

ZTT 

and replace the condition 7^ — 7^ G [a, [3) with the condition 5i + + 
■ ■ ■ + G [q!,/9) for 1 < i < M, A; > 0. The main difference is the 
absence of the 1/e in the logarithm. This is done so as to maintain a 
mean spacing tightly asymptotic to one. Set 

0(T)= 5^(7m-7.), 
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and sum by parts 



71 



Now, C{t) telescopes, and von Mangoldt's formula (j5(J|) implies that 
C{t) = t + 0(1), so that the r.h.s above equals N{T) + 0(log(T)), and 
6i is on average equal to one. In carrying out numerical experiments 
with zeros one can either use the normalization given in or 
For the theoretical purpose of examining leadings asymptotics of, say, 
the pair correlation, the factors appearing in these normalizations in 
the logarithm, l/(27re) or l/(27r), are not important as they only affect 
lower order terms. However, for the purpose of comparing numerical 
data to theoretical predictions it is crucial to include them. 

On a visit by Montgomery to the the Institute for Advanced Study, 
Freeman Dyson out that large unitary matrices have the same pair 
correlation. Let 

be the eigenvalues of a matrix in U(A^), sorted so that 

< 01 < ^2 • • • < ^iv < 27r. 
Normalize the eigenangles 
(56) = eiN/{2n) 

so that 6'j+i — 9i equals one on average. Then, a classic result in random 
matrix theory jM] asserts that 

N-'\{l<t<j<N,e,-eie[a,(3)}\ 

equals, when averaged according to Haar measure over U(A^) and let- 
ting N —>■ oo, 




( sin vrt ^ ^ 



dt. 



\ TXt 

Odlyzko [0] |02j has carried out numerics to verify Montgomery's con- 
jecture (jSni)- His most extensive data to date involves billions of zeros 
near the lO^^rd zero of C,i^s). With kind permission we reproduce 
Odlyzko's pair correlation picture in figure ^ 

This picture compares the l.h.s. of for many bins [a, 6) of size 
6 — a = .01 to the curve 

/ sinTrt^ ^ 

Odlyzko's histogram fits the theoretical prediction beautifully. Bogo- 
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Figure 1. The first graph depicts Odlyzko's pair cor- 
relation picture for 2 x 10^ zeros of ({s) near the lO^^rd 
zero. The second graph shows the difference between the 
histogram in the first graph and 1 — {{simrt) / (nt))'^ . In 
the interval displayed, the two agree to within about .002. 



molny and Keating |K] |BoKj . using conjectures of Hardy and Little- 



wood, have explained the role played by secondary terms in the pair 
correlation of the zeros of ({s) and these terms are related to ({s) on 
the one line. A nice description of these results are contained in |BK2j . 
Recently, Conrey and Snaith |CSj obtained the main and lower terms 
of the pair correlation using a conjecture for the full asymptotics of the 
average value of a ratio of four zeta functions rather than the Hardy- 
Littlewood conjectures. 

Montgomery's pair correlation theorem has been generalized 
by Rudnick and Sarnak (.RudSj to any primitive L-function, i.e. one 
which does not factor as a product of other L-functions, as well as to 
higher correlations which are defined in a way similar to the pair cor- 
relation. Again, there are severe restrictions on the fourier transform 
of the allowable test functions, and further, for L-functions of degree 
greater than three, Rudnick and Sarnak assume a weak form of the 
the Ramanujan conjectures. Bogomolny and Keating provide a heuris- 
tic derivation of the higher correlations of the zeros of ({s) using the 
Hardy-Littlewood conjectures |BoK2j . 

The author has tested the pair correlation conjecture for a number 
of L-functions. Figure El depicts the same experiment as in Odlyzko's 
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figure, but for various Diriclilet L-functions and L-functions associated 
to cusp forms. Altogether tfiere are eigliteen graphs. 

The first twelve graphs depict the pair correlation for all primitive 
Dirichlet L-functions, L{s, x) for conductors g = 3, 4, 5, 7, 8, 9, 11, 12, 
13, 15, 16, 17. Each graph shows the average pair correlation for each 
q, i.e. the pair correlation was computed individually for each x), 
and then averaged over x mod q. 

In the case of g = 3, 4 there is only one primitive L-f unction for either 
q, and approximately five million zeros were used for each (4, 772, 120 
and 5,003,411 zeros respectively to be precise). In the case of g = 
5, 7, 8, 9, 11, 12, 13, 15, 16, 17 there are 3, 5, 2, 4, 9, 1, 11, 3, 4, 15 primitive 
L-functions respectively. For g = 5, 7, 8, 9, 11, 12 either 2, 000, 000 ze- 
ros or 1,000,000 zeros were computed for each L{s,x), depending on 
whether x was real or complex. In the case of g = 16, 17 half as many 
zeros were computed. 

The last six graphs are for L-functions associated to cusp forms. The 
first of these six shows the pair correlation of the first 284, 410 zeros of 
the Ramanujan r L-function, corresponding to the cusp form of level 
one and weight twelve. The next five depict the pair correlation of the 
first 100, 000 zeros of the L-functions associated to the elliptic curves 
of conductors 11, 14, 15, 17, 19. These last six graphs use larger bins 
since data in these cases is more limited. 

The quality of the fit is comparable to what one finds with zeros of 
up to the same height. See, for example, figures 1 and 3 in |03j . It 
would be possible to extend the L(s, x) computations and obtain data 
near the lO^^th or higher zero, at least for reasonably sized g. Using 
the methods of section|21the time required to compute L{l/2 + it, x) is 
0{\qt\^/^), compared to 0(|t|^/2) for ((1/2 -Fit). Adapting the Odlyzko- 
Schonhage algorithm would allow for many evaluations of these L- 
functions at essentially the cost of a single evaluation. While such 
a computation might be manageable for Dirichlet L-functions, it is 
hopeless for cusp form L-functions where the time and also the number 
of Dirichlet coefficients required is 0{\N^^'^t\), i.e. linear in t. Here 
is the conductor of the L-function. Using present algorithms and 
hardware, it might be possible to extend these cusp form computations 
to t = 10^ or 10^. 

Slight care is needed to normalize these zeros correctly as the formula 
for the number of zeros of L{s) depends on the degree of the L-function 
and on its conductor. For Dirichlet L-functions L{s,x), X mod g, we 




Figure 2. Pair correlation for zeros of all primitive 
L{s,x), 3 < g < 17, the Ramanujan r L-function, and 
five elliptic curve L-functions 
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should normalize its zeros 1/2 + ^7 as follows: 

. log(|7|g/(27re)) 

^ = ^ 2^ 

For a cusp form L-function of conductor A^, we should take the following 
normalization: 

. log(|7|iV^/V(27re)) 
7 = 7 

TT 

From a graphical point of view, it is hard to display information 
concerning higher order correlations. Instead one can look at a statistic 
that involves knowing |KSj all the n-level correlations for characteristic 
functions, namely the nearest neighbour spacings distribution. 

In Figure El we display Odlyzko's picture for the distribution of the 
normalized spacings 6j for 2 x 10^ zeros of ({s) near the lO^^rd zero. 
This is computed by breaking up the into small bins and count- 

ing how many (5j's fall into each bin, and then comparing this against 
the nearest neighbour spacings distribution of the normalized eigenan- 
gles of matrices in U(iV), as — ^ oo, again averaged according to 
Haar measure on U(A^). The density function for this distribution is 
given |M] as 

n 

where A„(t) are the eigenfunctions of the integral operator 

(57) Mt)m = /' ""'7'"',^" /(y)^i/, 

J-1 T^{x-y) 

sorted according to 1 > Ao(t) > Ai(t) > ... > 0. See |U3j for a 
description of how the density function can be computed. 

In Figure|3]we display the nearest neighbour spacings distribution for 
the sets of zeros described above, namely millions of zeros of primitive 
L(s, x), with conductors 3 < g < 17, and hundreds of thousands of 
zeros of six cusp form L-functions. We also depict the nearest neighbour 
spacings for the first 500, 000 zeros of each of the 16 primitive x) 
with X mod 19 complex, and 1, 000, 000 zeros for the one primitive 
real x mod 19. 

Eight graphs are displayed. The first is for the 4, 772, 120 zeros of 
L{s.iX)^ X mod 3. The second one depicts the average spacings dis- 
tribution for all 76 primitive L(s, x), X mod q with 3 < g < 19, i.e. 
the spacings distribution was computed individually for each of these 
L-functions and then averaged. The next six graphs show the spac- 
ings distribution for the Ramanujan r L-function, and the L-functions 
associated to the elliptic curves of conductors 11, 14, 15, 17, 19. Again, 
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Figure 3. The first graph shows Odlyzko's nearest 
neighbour spacings distribution for 2 x 10^ zeros of ({s) 
near the lO^^rd zero. The second graph shows the differ- 
ence he computed between the histogram and the pre- 
dicted density function. Recently, Bogomolny, Bohigas 
and Leboeuf have explained the role of secondary terms 
in shaping the difference displayed. 



the fit is comparable to the fit one gets with the same number of zeros 
ofC(s). 

4.3. Density of zeros. Rather than look at statistics of a single L- 
function, we can form statistics involving a collection of L-functions. 
This has the advantage of allowing us to study the behaviour of our 
collection near the critical point where specific information about the 
collection may be revealed. This idea was formulated by Katz and 
Sarnak |KSj |KS2j who studied function field zeta functions and con- 
jectured that the various classical compact groups should be relevant 
to questions about L-functions. 

While the eigenvalues of matrices in all the classical compact groups 
share, on average, the same limiting correlations and spacings distri- 
butions, their characteristic polynomials do exhibit distinct behaviour 
near the point z = 1. Using the idea that the unit circle for char- 
acteristic polynomials in the classical compact groups correponds to 
the critical line, with the point 2; = 1 on the unit circle corresponding 
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Figure 4. Nearest neighbour spacings distribution for 
several Dirichlet and cusp form L-f unctions. The first 
is for L(s, x), q — 3. The second is the average nearest 
neighbour spacing for all primitive L{s,x), 3 < g < 19. 
The last six arc for the Ramanujan r L- function, and 
five L-functions associated to elliptic curves. 
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to the critical point, Katz and Sarnak were led to formulate conjec- 
tures regarding the density of zeros near the critical point for various 
collections of L-functions. This is detailed in section 14.3.11 below. 

The fact that different families of L-funtions exhibit distinct be- 
haviour near the critical point is illustrated in figure This plot 
depicts the imaginary parts of the zeros of many L{s, x) with x ^ 
generic non-real primitive Dirichlet character for the modulus q, with 
5 < q < 10000. Other than the fact that, at a fixed height, the ze- 
ros become more dense proportionally to log q, the zeros appear to be 
uniformly dense. 

This contrasts sharply with the plot in figure IHl which depicts the 
zeros of L{s, Xd) where Xd is a real primitive character (the Kronecker 
symbol), and d ranges over fundamental discriminants with —20000 < 
d < 20000. Here we see the density of zeros fluctuating as one moves 
away from the real axis. 

Other features can be seen in the plot. First, from the white band 
near the x-axis we notice that the lowest zero for each L{s, Xd) tends 
to stay away from the critical point. We can also see the effect of 
secondary terms on this repulsion. The lowest zero for d > tends to 
be higher than the lowest zero for d < 0. This turns out to be related 
to the fact that the F-factor in the functional equation for L{s, Xd) is 
F(s/2) if > 0, but is F((s + l)/2) when d < 0. 

We can also see slightly darker regions appearing in horizontal strips. 
The first one occurs roughly at height 7., half the height of the first zero 
of C{s)- These horizontal strips are due to secondary terms in the den- 
sity of zeros for this collection of L-functions which include |Ru3j |CSj 
a term that is proportional to 

C(l + 2zt) ■ 

This is large when C(l + 2«t) is small. Surprisingly, (^(1 + iy) and 
C,{\/2 + iy) track each other very closely, see figure [71 and the minima 
oi \C,{1 + iy)\ appear close to the zeros oi (,{1/2 + iy). This is similar 
to a phenomenon that occurs when we look at secondary terms in the 
pair correlation of the zero of which also involves C(s) on the one 
line |BK2] . 

4.3.1. n-level density. The n-level density is used to measure the av- 
erage density of the zeros of a family of L-functions or matrices. It is 
arranged to be sensitive to the low lying zeros in the family, i.e. those 
near the critical point if we are dealing with L-functions, and those near 
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Figure 5. Zeros of L(s,x) with x a generic non-real 
primitive Dirichlet character for the modulus q, with 5 < 
q < 10000. The horizontal axis is q and, for each L[s, x), 
the imaginary parts of its zeros up to height 15 are listed. 



the point 2; = 1 on the unit circle if we are dealing with characteristic 
polynomials from the classical compact groups. 

Let A be an X matrix in on of the classical compact groups. 
Write the eigenvalues of A as Xj = e*^J with 

< ^1 < . . . < ^jv < 27r. 

Let 

H^-\A, f)= J2 f (^nN/{2n), 9,„N/ {2n)) 

l<jl,..-,in<JV 
distinct 

with / : M" — > M, bounded, Borel measurable, and compactly sup- 
ported. Because of the normalization by N/{2tt), and the assumption 
that / has compact support, H^"'\A, f) only depends on the small 9/s. 
Katz and Sarnak KS proved the following family dependent result: 



/• /'OO /"CXD 

(58) lim / H^''\AJ)dA= / ... / W^\x)f{x)dx 

Jg(n) Jo Jo 



'G(iV) 

for the following families: 
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Figure 6. Zeros of L{s,Xd) with Xd(n) = (^), the Kro- 
necker symboL We restrict d to fundamental discrimi- 
nants -20000 < d < 20000. The horizontal axis is d 
and, for each L{s,Xd), the imaginary parts of its zeros 
up to height 30 are listed. A higher resolution image can 
be obtained from the author's webpage under 'Publica- 
tions'. 
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Figure 7. A graph illustrating that, at least initially, 
the minima of + iy)\ occur very close the zeros of 
|C(l/2 + iy)\- The dashed line is the graph of the former 
while the solid line is the graph of the the latter. 



G 




U(iV),U.(iV) 


det {Kf){xj,Xk))i<j<n 

l<k<7i 


USp(A^) 


det {K_i{xj,Xk))i<]<n 

l<k<n 


S0(2A^) 


det {Ki{xj,Xk))i<j<n 

l<k<n 


S0(2A^+ 1) 


det {K_i{xj,Xk)) i<j<n + Zl"=i ^i^i^) 'i^t {K^i{xj, Xk)) i<,^.<n 


with 



sin(7r(a; — y)) sin(7r(a; + |/)) 
7r(x — y) Ti{x + y) 



Here 

U^(A^) = {Ae \J{N) : det(A)" = 1}. 
The delta functions in the S0(2A^ + 1) case are accounted for by the 
eigenvalue at 1. Removing this zero from ()58j) yields the same Wq"^ as 
for USp. 
Let 

D{X) = {d a, fundamental discriminant : \d\ < X} 

and let Xdl^) ~ (n) Kronecker's symbol. Write the non-trivial zeros 
of L(s, Xd) as 

l/2 + i^f\ j = ±l,±2,... 
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sorted by increasing imaginary part, and 




The author proved |R,u2j that 




) 



(59) 



poo poo 

/ . . . / f{x)wiS^{x)dx, 
Jo Jo 



where 




supported in Yl^=i I'^^^l < 1- This generahzed the n = 1 case that had 
been achieved earher |()zSj |KS2j . Assuming the Riemann Hypothesis 
for all L{s,Xd), the n = 1 case has been extended to / supported in 
(—2,2) [OzS2j [KS3 . Chris Hughes has an alternate derivation of (j^IJj) 
appearing in the notes of these proceedings. 

This result confirms the connection between zeros of L{s, Xd) and 
eigenvalues of unitary symplectic matrices and explains the repulsion 
away from the critical point and the fluctuations seen in figure IHl at 
least near the real axis, because, when n = 1, the density of zeros is 
described by the function wjj^ (x) which equals 



At height x, we therefore also expect, as we average over larger and 
larger \d\, for the fluctuations to diminish proportional to 1/x. How- 
ever, if we allow x to grow with d then the fluctuations actually persist 
due to secondary fluctuating terms that can be large if x is allowed to 
grow with d |Ru3j |CSj . 

The above suggests that the distribution of the lowest zero, i.e. the 
one with smallest imaginary part, in this family of L-functions ought 
to be modeled by the distribution of the smallest eigenangle of char- 
acteristic polynomials in USp(A^), with — > oo. Similary we expect 
that the distribution, say, of the second lowest zero ought to fit the 
distribution of the second smallest eigenangle. 

The probability densities describing the distribution of the smallest 
and second smallest eigenangles, normalized by N/{2tt), for charac- 
teristic polynomials in USp(A^), with A^ even and tending to oo are 



1 - 



sin(27ra;) 



271X 
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given |KSj respectively by 



and 
where 



z.i(USp)(t) = -^i?_,o(t) 
i.2(USp)(t) = -|(E_,o(t) + i?-,i(t)), 



^-,o(t) = n(l-^2,+i(2t)) 



j=0 



k=l J=o 



Here, the Aj(t)'s are the eigenvalues of the integral equation in (jFTjl . 

This also suggests that the means of the the first and second lowest 
zeros are given by 

1 Z""" 
i^oofmY^ E ^^'^^^= / t^2(USp)(t)cit = 1.76... 
' ^ deD{x) •'^ 

However, the convergence to the predicted means is logarithmically 
slow due to secondary terms of size 0(1/ log(X)). Consequently, when 
comparing against the random matrix theory predictions, one gets a 
better fit by making sure the lowest zero has the correct mean. This 
can be achieved by rescaling the data, further multiplying, for a set D 
of fundamental discriminants, '^^i^'ld by 



(60) ■78(^^E^f''^^ 

and '-^fhd by 



deD 



(61) 1-76 r^E^?^^" 

V ' d&D 

In figures IHl and ini we use the normalization described above. For our 
data set, the denominator in (p?j) equals .83, and, in equals 1.84. 

In figure IHl we depict the 1-level density of the zeros of L(s,Xd) for 
7243 prime \d\ lying in the interval (10^^ 10^^ _^ 200000). These zeros 
were computed in 1996 as part of the authors PhD thesis |Ruj . Here we 
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divide the x-axis into small bins, count how many normalized zeros of 
L{s, Xd) lie in each bin, divide that count by the number of namely 
7243, and compare that to the graph of 1 — sin(27rx)/ {2txx). 

1.5 -| 




0.0 0.5 1.0 1.5 2.0 



Figure 8. Density of zeros of L{s,Xd) for 7243 prime 
values of \d\ lying in the interval (10^^ 10^^ ^ 200000). 
Compared against the random matrix theory prediction, 
1 — sin(27ra;)/(27ra;). 

In figure ini we depict the distribution of the lowest and second lowest 
normalized zero for the set of zeros just described. These are com- 
pared against vi and U2 which were computed using the same program, 
obtained from Andrew Odlyzko, that was used in [03 . 

In figure ^1 we depict the 1-level density and distribution of the low- 
est zeros for quadratic twists of the Ramanujan r L-function, Lt-{s, Xd), 
d> Q. For this family of L-functions, one can prove |Ru2j a result sim- 
ilar to (jSni) but with W^usp replaced with M/so(even), and the support 
of / reduced to X]r=i 1"^*! ^ -'-Z^- '^^^ 1-level density is therefore given 
by 1 + sin(27rx)/(27rx) and the probability density for the distribu- 
tion of the smallest eigenangle, normalized by 2N/ (27r), for matrices in 
S0(2iV), with ^ oo, is given [KB] by 

, oo 

i.i(SO(even))(t) = --n(l-A2,(2t)), 

j=0 

whose mean is .32. the figure uses 11464 prime values of \d\ lying 
in (350000,650000), and the zeros were normalized by 21^-, and then 
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Figure 9. Distribution of the lowest and second low- 
est zero of L{s,Xd) for 7243 prime values oi \d\ lying in 
the interval (10^^ 10^^ _^ 200000). Compared against the 
random matrix theory predictions. 

rescaled so as to have mean .32 rather than .29. The choice of using 
2ld for normalizing the zeros is the correct one up to leading term, 
but is slightly adhoc and by now a better understanding of a tighter 
normalization up to lower terms has emerged |CFKRS^ |Ru3j. 

4.4. Value distribution of L-functions. Keating and Snaith initi- 
ated the use of random matrix theory to study the value distribution 
of L-functions with their important paper |KeSj where they consider 
moments of characteristic polynomials of unitary matrices and con- 
jecture the leading-order asymptotics for the moments of C,{s) on the 
critical line. This was followed by a second paper |KeS2j along with 
a paper by Conrey and Farmer jCF| which provide conjectures for the 
leading-order asymptotics of moments of various families of L-functions 
by examining analogous questions for characteristic polynomials of the 
various classical compact groups. 

Keating and Snaith's technically impressive work also represents a 
philosophical breakthrough. Until their paper appeared, one would 
compare, say, statistics involving zeros of C,{s) to similar statistics for 
eigenvalues of N x N unitary matrices, with — ^ oo. However, their 
work compares the average value of ({1/2 + it) to the average value of 
N X N unitary characteristic polynomials evaluated on the unit circle, 
with N ~ log(t/(27r)). This choice of is motivated by comparing 
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Figure 10. One-level density and distribution of the 
lowest zero of even quadratic twists of the Ramanujan 
T L-function, LT-{s,Xd), for 11464 prime values of d > 
lying in the interval (350000,650000). 



local spacings of zeros, for example v.s. A slightly differ- 

ent approach to this choice of proceeds by comparing functional 
equations of L-functions to functional equations of characteristic poly- 
nomials jOEKES]- 

At first sight, it seems strange to compare the Riemann zeta function 
which has infinitely many zeros to characteristic polynomials of finite 
size matrices. However, this suggests that a given height t, the Riemann 
zeta function can be modeled locally by just a small number of zeros, 
as well as by more global information that incorporates the role played 
by primes. Recently, Gonek, Hughes, and Keating have developed such 
a model [UHKj. 

Below we describe three three specific examples where random ma- 
trix theory has led to important advances in our understanding of the 
value distribution of L-functions. These concern the families: 

(1) ({1/2 + it), where we average over t. 

(2) L(l/2, Xd), where we average over fundamental discriminants d. 

(3) L^(l/2,Xd), quadratic twists of the L-function associated to 
an elliptic curve E over Q, where we average over fundamental 
discriminants d. 
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These three are examples of unitary, unitary symplectic, and even or- 
thogonal families respectively |KS2j |CFKRSj . Note that in the last 
example, we normalize the Dirichlet coefficient of the L-function as 
in (jHTj) so that the functional equation of L£;(l/2, x^) brings s into 
1 — s with the critical point being s = 1/2. 

We first illustrate that these three examples exhibit distinct be- 
haviour by contrasting their value distributions. The first-order asymp- 
totics for the moments of |C(l/2 + it)\ are conjectured by Keating and 
Snaith |KeSj to be given by 

(62) ^ IC(l/2 + tt) rdt ^ a^i, n r(j + 1/2)'^^ ' J?r>-1, 

with ~ log(T) and 0^/2 defined by (|T!?jl . 

For quadratic Dirichlet L-functions Keating and Snaith |KeS2j con- 
jecture that 
(63) 

1 V- y . .2iv.A r(iv + j + i)r(j + i + r) 

^ ' ^ ' r(iv + . + i + r)r(, + i)' > 

with ~ log(X)/2, where the sum runs over fundamental discrimi- 
nants \d\ < X, and, as suggested by Conrey and Farmer |CFj . 



V 



Next, let q be the conductor of the elliptic curve E. Averaging 
over fundamental discriminants and restricting to discriminants for 
which LE{s,Xd) has an even functional equation, the conjecture as- 
serts [CF] |K^ that 
(64) 

V LMi/2, - c.2^^- n + ' - '^'^^\-^^ + , ^r>-i/2, 



I^WI ^^'"^"^ |ir(Ar + ,-i + r)r(j-i; 

(d,g) = l 
even funct cqn 

with A^ ~ log(X) and 

/ ^Nfc(fc-l)/2 
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In the above equation, ap stands for the pth coefficient of the Dirichlet 
series of Le- 

In the case of the Riemann zeta function we take absolute values, 
|C(l/2 + it)\, otherwise the moments would be zero. In the other two 
cases, the L-values are conjectured to be non-negative real numbers, 
hence we directly take their moments. 

We should also observe that, while statistics such as the pair correla- 
tion or density of zeros involving zeros of L-functions have arithmetic 
information appearing in the secondary terms, moments already re- 
veal such behaviour at the level of the main term. This reflects the 
global nature of the moment statistic as compared to the local nature 
of statistics of zeros that have been discussed. 

Using the above conjectured asymptotics we can naively plot value 
distributions. Figure ^2 compares numerical value distributions for 
data in these three examples against the counterpart densities from 
random matrix theory. Notice that these three graphs behave distinctly 
near the origin. The solid curves are computed by taking inverse Mellin 
transforms, as in (fT^ . of the right hand sides of equations (jHSl), 
and but without the arithmetic factors a^, bk, Ck- Shifting the 
inverse Mellin transform line integral to the left, the location of the 
first pole in each integrand dictates the behaviour of the corresponding 
density functions near the origin. The locations of these three poles are 
at r = — 1, —3/2, and —1/2 respectively. Taking t to be the horizontal 
axis, near the origin the first density is proportional to a constant, the 
second to t^^^, and the third to In forming these graphs one takes 

as described above, so that the proportionality constants do depend 
on N. As N grows, these graphs tend to get flatter. 

The first graph is reproduced from jKeSj . In the second and third 
graphs displayed, a slight cheat was used to get a better fit. The 
histograms were rescaled linearly along both axis until the histogram 
matched up nicely with the solid curves. We must ignore the arith- 
metic factors when taking inverse Mellin transforms since these factors 
are known |CGoj to be functions of order two and cause the inverse 



Mellin transforms to diverge. To properly plot the correct value distri- 
butions we would need to use more than just the leading-order asymp- 
totics. Presently, our knowledge of the moments of various families of 
L-functions extends beyond the first-order asymptotics, but only for 
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positive integer values of r (even integer in the case of |C(l/2 + it)\'^), 
however, one would need to apply full asymptotics for complex val- 
ues of r. The paper by Conrey, Farmer, Keating, the author, and 
Snaith |CFKRSj conjectures the full asymptotics, for example, of the 
three moment problems above, but for integer r, with corresponding 
theorems in random matrix theory given in [CFKRS2]. The paper of 
Conrey, Farmer and Zirnbauer goes even beyond this stating conjec- 
tures for the full asymptotics of moments of ratios of L-functions, and, 
using methods from supersymmetry, proving corresponding theorems 
in random matrix theory |CFZj . Another paper, by Conrey, Forrester, 
and Snaith, uses orthognal polynomials to obtain alternative proofs of 
the random matrix theory theorems for ratios |CFSj . 



4.4.1. Moments of\({l/2 + it)\. Next we describe the full moment con- 
jecture from ICFKR^Sj for |C(l/2 + it)\. In that paper, the conjecture 
is derived heuristically by looking at products of zetas shifted slightly 
away from the critical line and then setting the shifts equal to zero. 

The formula is written in terms of contour integrals and involve the 
Vandermonde: 



A(zi, ...,Zm)= Y\. (^i " ^»)- 

l<i<j<m 



Suppose g(t) = f(t/T) with / : M"*" —>■ M non- negative, bounded, 
and integrable. The conjecture of |CFKRSj states that, as T cxo, 

/»00 /"OO 

/ \ai/2 + tt)\'' git) dtr^ / P,i\ogit/{2n)))git)dt, 
Jo Jo 

where is the polynomial of degree fc^ given by the 2fc-fold residue 



Pk(x) 



(27ri) 



2k 



2k 

n 



dzi... dz 



2k, 



,2k 



where one integrates over small circles about Zi = 0, with 



k k 

G{zi, . . .,Z2k) = Ak{Zi, . . • ,^2fc) ]^]^C(1 + Zi- Zk+j), 

i=l j=l 
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Figure 1 1 . Value distribution of L-functions compared 
to the random matrix theory counterparts. The first pic- 
ture, depicts the value distribution of |^(l/2 + it)|, with 
t near 10^ the second of L(l/2, Xd) with 800000 < \d\ < 
10^ and the third of LE,^{l/2,Xd), with -85000000 < 
d<0, d= 2,6,7,8, 10 mod 11. 



and is the Euler product 

k k , ^ \ /"l / 

Aw=nnn(i-3T;^)/n(i 

p i=i 0=1 V P ' J Jo ,=1 V 



P 



2^ ^3 



-1 



p2 



p 



n n I _ p^fc+i 

p m=l i^m 



de 
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When A; = 1 or 2, this conjecture agrees with theorems for the 
full asymptotics as worked out by Ingham P and Heath-Brown re- 
spectively jHj. In the first case Ai{z) = 1 and in the second case 
^2(-2) = C(2 + Zi + Z2 — Z3 — Z4)~^, and one can write down the coef- 
ficients of the polynomials Pk{x) in terms of known constants. When 
k = 3 the product over primes becomes rather complicated. However, 
one can numerically evaluate |CFKRS3] the coefficients of Psix) and 



the polynomial is given by: 

Psix) = 0.000005708527034652788398376841445252313 

+0.00040502133088411440331215332025984x^ 
+0.011072455215246998350410400826667x^ 
+0.14840073080150272680851401518774x^ 
+1.0459251779054883439385323798059x5 
+3.984385094823534724747964073429 
+8.60731914578120675614834763629 x^ 
+10.274330830703446134183009522x2 
+6.59391302064975810465713392 x 
+0.9165155076378930590178543. 

In the k = 3 case the moments of |C(l/2 + it)\ have not been proven, 
and it makes sense to test the moment conjecture numerically. Table El 
reproduced from ^CFKRS], depicts 

(65) r \ai/2 + tt)\'dt 



as compared to 



c 



(66) / P;{\og{t/27r))dt, 

Jc 

along with their ratio, for various blocks [C, D] of length 50000, as well 
as a larger block of length 2,350,000. 

4.4.2. Moments of L{l/2,Xd)- Another conjecture listed in |CFKRSj 
concerns the full asymptotics for the moments of L(l/2, Xd)- We quote 
the conjecture here: 

Suppose g{t) = f(t/T) with / : IR+ M non- negative, bounded, 
and integrable. Let Xa{s) = \d\2~^X{s,a) where a = if c? > and 
a = 1 if d < 0, and 

X{s,a) = 7r^^-^T{ ] /T ' 
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[C, D] 


conjecture (j6SI) 


reality (jSSj) 


ratio 


[0,50000] 


7236872972.7 


7231005642.3 


.999189 


[50000,100000] 


15696470555 


3 


15723919113 


6 


1.001749 


[100000,150000] 


21568672884 


1 


21536840937 


9 


.998524 


[150000,200000] 


26381397608 


2 


26246250354 


1 


.994877 


[200000,250000] 


30556177136 


5 


30692229217 


8 


1.004453 


[250000,300000] 


34290291841 





34414329738 


9 


1.003617 


[300000,350000] 


37695829854 


3 


37683495193 





.999673 


[350000,400000] 


40843941365 


7 


40566252008 


5 


.993201 


[400000,450000] 


43783216365 


2 


43907511751 


1 


1.002839 


[450000,500000] 


46548617846 


7 


46531247056 


9 


.999627 


[500000,550000] 


49166313161 


9 


49136264678 


2 


.999389 


[550000,600000] 


51656498739 


2 


51744796875 





1.001709 


[600000,650000] 


54035153255 


1 


53962410634 


2 


.998654 


[650000,700000] 


56315178564 


8 


56541799179 


3 


1.004024 


[700000,750000] 


58507171421 


6 


58365383245 


2 


.997577 


[750000,800000] 


60619962488 


2 


60870809317 


1 


1.004138 


[800000,850000] 


62661003164 


6 


62765220708 


6 


1.001663 


[850000,900000] 


64636649728 





64227164326 


1 


.993665 


[900000,950000] 


66552376294 


2 


65994874052 


2 


.991623 


[950000,1000000] 


68412937271 


4 


68961125079 


8 


1.008013 


[1000000,1050000] 


70222493232 


7 


70233393177 





1.000155 


[1050000,1100000] 


71984709805 


4 


72919426905 


7 


1.012985 


[1100000,1150000] 


73702836332 


4 


72567024812 


4 


.984589 


[1150000,1200000] 


75379769148 


4 


76267763314 


7 


1.011780 


[1200000,1250000] 


77018102997 


5 


76750297112 


6 


.996523 


[1250000,1300000] 


78620173202 


6 


78315210623 


9 


.996121 


[1300000,1350000] 


80188090542 


5 


80320710380 


9 


1.001654 


[1350000,1400000] 


81723770322 


2 


80767881132 


6 


.988303 


[1400000,1450000] 


83228956776 


3 


83782957374 


3 


1.006656 


[0,2350000] 


3317437762612.4 


3317496016044.9 


1.000017 



Table 3. Sixth moment of C, versus the conjecture. The 
'reahty' column, i.e. integrals involving (, were com- 
puted using Mathematica. 



That is, Xd{s) is the factor in the functional equation L{s, Xd) = 
Xd{s)L{l — s,Xd)- Summing over negative fundamental discriminants 
d we have, as T ^ oo. 



d<0 d<0 
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where Qk is the polynomial of degree k{k + l)/2 given by the /c-fold 
residue 

G'(^l,...,^fc)A(^l^■■■,^^)^ j^fe 

~k 

where 

k 

i=i i<i<i<fc 

and is the Euler product, absolutely convergent for |3?2j| < 
defined by 

5.(^1,...,^.) =n n (i-^iT^) 

x(5(nO--W) '^n(i-^-W) 

We can also sum over d > but then need to replace X(l/2 + Zj^ 1) 
with X(l/2 + Zj,0). 

This conjecture agrees with theorems in the case of A; = 1, 2, 3 [J] [S] 
(only the leading term has been checked in the case of /c = 3, but in 
principle the lower terms could be verified). 

Figure [T21 reproduced from |CFKRSj . depicts, for /c = 1, . . . , 8 and 
X = 10000, 20000,..., 10^ 

0<d<X 

divided by 

0<d<X 

4.4.3. Vanishing of LE{l/2,Xd)- In |(^KRSj . Conrey, Keating, the au- 
thor, and Snaith apply the moment conjecture (jMj) to the problem 
of predicting asymptotically the number of vanishings of LE{l/2,Xd)- 
Using the fact that these L-values are discretized, for example via 
the Birch and Swinnerton-Dyer conjecture or the theorem of Kohnen- 
Zagier |KZj , and by studying, up to leading term and for small values. 



Qk{x) 



k\ (2^ 



68 



MICHAEL RUBINSTEIN 




Figure 12. Horizontal axis in each graph is X. These 
graphs depict the first eight moments, sharp cutoff, of 
L(l/2,Xrf), < d < X divided by the conjectured 
value, sampled at X = 10000, 20000, . . . , 10^. We see 
the graphs fluctuating above and below one. Notice that 
the vertical scale varies from graph to graph 
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the density function predicted by they conjectured that 

J2 l~«£;X3/4log(Xf^. 

deD(X) 

(d,9) = l 
even funct cqn 
i'i5(l/2,Xd)=0 

The power on the logarithm depends on the underlying curve E be- 
cause, in the Birch Swinnerton Dyer conjecture, the Tamagawa factors 
can contribute powers of 2 depending on the prime factors of d and 
on E and this affects the discretization. The constant as depends on 
c_i/2 and the real period of E, but also on some extra subtle arithmetic 
information that seems to be related to Delaunay's heuristics for Tate- 
Shafarevich groups [DeJ and is not yet fully understood. Numerical 
evidence in favour of this conjecture is presented in |CKRS2j . One can 
skirt these delicate issues, the power on the logarithm and the constant 
ttE, as follows. 

Let p \ qhe prime. Sort the d's for which L£;(l/2, Xd) = by residue 
classes mod p, according to whether Xd{p) = 1 or —1, and consider the 
ratio 

dGD(X) 1 

(d,9) = l 
even funct eqn 
i'B(l/2,Xd)=0 



Rp{X) 



Xd(p)=i 



(d,g) = l 
even funct eqn 
iB{l/2.Xd)=0 

Xd(p)=-i 



One can formulate |CKRS2j |CFKRSj conjectures for the moments in 



these two subfamilies and the moments agree except for a factor that 
depends on p. By considering this ratio, the powers of X, of logX, and 
the constant should all cancel out, except for a single factor that 
depends on p. This leads to a conjecture ^CKRSj for Rp{X): 



Rp = lim Rp{X) 



IP + 1 — a. 



v 

5 



x-^oo -^^ y p + 1 + a 

where Op denotes the pth coefficient of the Dirichlet series for Le- The 
square root in this conjecture is a consequence of the moments having 
a pole at r = — 1/2. 

We end this paper with a plot that substantiates this conjecture. 
Figure [121 compares, for one hundred elliptic curves E, the predicted 
value of Rp to the actual value Rp{X), with X = 10^ and the set of 
(i's restricted to certain residue classes depending on E as described 
in )CKRS2j . The L- values were computed in this special case by ex- 
ploiting their connection to the coefficients of certain weight three 
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halves modular forms and using a table of Rodriguez- Villegas and 
Tornaria |Rl'j . 

The horizontal axis is p. For each p, and each of the one hundred el- 
liptic curves E we plot Rp{X) — Rp. We see the values fluctuating about 
zero, most of the time agreeing to within about two percent. The con- 
vergence in X is predicted from secondary terms to be logarithmically 
slow and one gets a better fit by including more terms jCKRS2j . 
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Figure 13. A plot for one hundered elliptic curves of 
Rp{X) -Rpioi2<p< 3500, X = 10^ 
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